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Abstract 

First Evidence For Atmospheric Neutrino-Induced Cascades with the IceCube 

Detector 

by 

Michelangelo Vincent DAgostino 
Doctor of Philosophy in Physics 

University of California, Berkeley 
Professor P. Buford Price, Chair 

IceCube is an all-flavor, cubic kilometer neutrino telescope currently under construction 
in the deep glacial ice at the South Pole. Its embedded optical sensors detect Cherenkov 
light from charged particles produced in neutrino interactions in the ice. For several years 
IceCube has been detecting muon tracks from charged-current muon neutrino interactions. 
However, IceCube has yet to observe the electromagnetic or hadronic particle showers or 
"cascades" initiated by charged-current or neutral-current neutrino interactions. The first 
detection of such an event signature is expected to come from the known flux of atmospheric 
electron and muon neutrinos. 

A search for atmospheric neutrino-induced cascades was performed using 275.46 
days of data from IceCube's 22-string configuration. Reconstruction and background re- 
jection techniques were developed to reach, for the first time, a signal-to-background ratio 
~1. Above a reconstructed energy of 5 TeV, 12 candidate events were observed in the full 
dataset. The signal expectation from the canonical Bartol atmospheric neutrino flux model 
is 5.63 ± 2.25 events, while the expectation from the atmospheric neutrino flux as measured 
by IceCube's predecessor array AMANDA is 7.48 ± 1.50 events. Quoted errors include the 
uncertainty on the flux only. 

While a conclusive detection can not yet be claimed because of a lack of background 
Monte Carlo statistics, the evidence that we are at the level of background suppression 
needed to see atmospheric neutrino-induced cascades is strong. In addition, one extremely 
interesting candidate event of energy 133 TeV survives all cuts and shows an intriguing 



2 

double pulse structure in its waveforms that may signal the "double bang" of a tau neutrino 
interaction. 



Professor P. Buford Price 
Dissertation Committee Chair 
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Chapter 1 

Neutrinos in Physics and 
Astrophysics 

1.1 Introduction: Why Neutrinos? 

For the last half-century, scientists have been traveling to remote, often harsh 
locations to capture neutrinos. They've spent countless days in deep underground mines 
[1, 2, 3, 4] and caverns situated next to nuclear reactor cores [5, 6, 7]. They've ventured 
to frozen Siberian lakes [8] and to the Amundsen-Scott South Pole Station [9, 10] . They've 
even endured the inhospitable climate of the French Riviera [11]. So what's all the fuss 
about? 

Over that same time period, the neutrino has moved from a minor player, invented 
to save energy conservation but never expected to be detected, to a central research topic 
in particle physics and a hopeful, emerging research topic in astrophysics. The discovery 
of neutrino oscillations has opened the door onto a new and as-yet unknown mechanism 
beyond the Standard Model that must be responsible for giving neutrinos their tiny masses. 
Next-generation oscillation experiments hope to probe CP violation in the neutrino sector, 
which many think could explain the matter-antimatter asymmetry of the universe. In the 
realm of astrophysics, neutrinos have confirmed our understanding of solar physics and 
given us a small glimpse into supernova explosions. The scientists who are building large 
detectors to detect supernova neutrinos and high energy astrophysical neutrinos are hoping 
to shed light on some of the processes at work inside stellar explosions, supermassive black 
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holes, gamma ray bursts, and other violent astrophysical objects. 

This dissertation is concerned with a small piece of this exciting tapestry of neu- 
trino physics and astrophysics. The IceCube detector is a next-generation, all-flavor neu- 
trino telescope currently under construction in the deep glacial ice at the South Pole. Ice- 
Cube's primary goal is to observe astrophysical neutrinos and, from them, to discover the 
sources of the highest energy cosmic rays. In order to achieve these goals, IceCube must 
prove that it can observe the various interaction signatures of different neutrino flavors using 
the only calibration source of high energy neutrinos that we have — atmospheric neutrinos. 
In particular, this dissertation focuses on a search for neutrino-induced particle showers or 
"cascades" from atmospheric neutrinos. 

1.2 Atmospheric Neutrinos 

The earth is constantly bombarded by a stream of high energy particles from 
space — the so-called cosmic rays. While the composition of cosmic-rays is energy dependent, 
roughly 75% are protons, 15% are helium, and 10% are heavier nuclei like carbon, oxygen, 
and nitrogen up through iron [12]. Cosmic rays reach earth with a tremendous amount 
of energy. The spectrum is a power law, 4^ ~ E~ 2 ' 7 up to energies around 10 3 TeV 
(the "knee"), where the spectrum steepens to ^ ~ E~ 3 . At an energy around 10 6 TeV 
(the "ankle" ) the spectrum hardens again to ^ ~ E~ 2 ' 7 before suppression begins around 
4 x 10 7 TeV [13]. This suppression is due to interactions of the highest energy cosmic rays 
with the photons of the cosmic microwave background, the so-called GZK effect [14, 15]. 
The ankle is thought to signal a transition from galactic to extragalactic sources. Figure 1.1 
shows the measured cosmic ray energy spectrum. 

Interactions of these primary cosmic ray nuclei with the atmosphere produce un- 
stable mesons which can further interact or decay to secondary muons, muon neutrinos, and 
electron neutrinos. Since this dissertation is concerned with a search for neutrino-induced 
cascades from atmospheric neutrinos, we'll spend some time discussing a few salient features 
of the atmospheric neutrino spectrum. 

The atmospheric neutrino spectrum is made up of contributions from the decay 
of muons, pions, and kaons, and its makeup is determined by the competition between 
interaction and decay. To quantify this, we follow the atmospheric model of [12]. The 
pressure is proportional to vertical slant depth X v , where X v is measured from the top of 
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Figure 1.1: Measured cosmic ray energy spectrum. From [16]. 
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the atmosphere down in units of g/cm 2 . The density is p = —dX v /dh, where h is the height 
in m measured from the ground up. We can write 

P X v 

— = o£ / 

p —dX v /dh 

For an isothermal atmosphere, 

X v = X e- h/ho 

where Xq ~ 1030 g/cm 2 is the total atmospheric vertical slant depth. Both the temperature 
and the scale height actually decrease with altitude, but at sea level ho ~ 8.4 km. This 
model for the atmosphere implies that 

p = —dX v /dh = X v /ho w — 

ho 

where X is the slant depth along the particle's trajectory measured from the top of the 
atmosphere down and 8 is the particle's zenith angle. We want to find the particle decay 
length d in g/cm 2 . This is given by 



E XcosO EX cos 6 EX cos 9 
d = jcrp = — g ct — = — = 

mc z riQ mc z ho/cT e 

where the critical energy e = mc ci h ° = The critical energy is the energy at which the 
particle decay length is equal to the scale height of the atmosphere. The critical energies 
for particles important for atmospheric neutrinos are given in table 1.1. 

At low energies such that E <S e/cos#, the decay length becomes very short. 
Particles will decay before they have a chance to interact. The neutrinos from these decays 
will have an energy spectrum that follows the mesons and primary cosmic rays. At energies 
E S> e/cos#, the decay length becomes very long. Particles are likely to interact before 
they can decay. The neutrinos will have an energy spectrum that is one power steeper than 
that of the mesons and primary cosmic rays because of the decay length dependence on 
E. The l/cos# dependence also gives atmospheric neutrinos a characteristic zenith angle 
dependence which peaks towards the horizon, where decay lengths are the longest. 

The most important unstable mesons produced in cosmic ray interactions are the 
charged pions and the charged and neutral kaons. The decay of pions and kaons produces 
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Table 1.1: Critical energies of various particles. 



Particle 


e (GeV) 




1.0 




115 


K± 


850 


K° L 


205 


D ± 


4.3 x 10 7 


D° 


9.2 x 10 7 



the bulk of atmospheric neutrinos. The pion decay (essentially 100% branching ratio) [17] 
is 

I 

At low neutrino energies, this decay dominates the spectrum. We therefore expect two muon 
neutrinos for every electron neutrino at low energies. As energy increases, however, the 
muon becomes increasingly boosted and eventually lives long enough to hit the surface of the 
earth before decaying. At this point, we lose muon decay as a source of electron neutrinos. 
The transition happens first for vertically downgoing muons and last for horizontal muons, 
which have a longer path over which to decay. 

The corresponding charged kaon decay (63.4% branching ratio [17]) is 

K ± ->n ± + z^(z^) 

Since the critical energy for kaons is higher than that for pions (850 GeV vs. 115 GeV), 
the atmospheric neutrino spectrum becomes dominated by kaon decay at around 100 GeV. 
These are the energies most relevant for IceCube. 

The charged kaon decays don't produce an appreciable flux of electron neutrinos. 



G 



The main source of high energy electron neutrinos is the small contribution from the neutral 
kaon K(3 decay mode 

Kl n ± + + v e {v e ) 

This has a branching ratio of 40.5% [17]. Also, has a lower critical energy than the 
charged kaons (205 GeV vs. 850 GeV). At high energies, we therefore expect the electron 
neutrino flux to be about a factor of 20 below the muon neutrino flux [18] . 

The charged pion and charged kaon decays are two-body processes, so in the 
rest frame of the decaying parent the energy and momentum of the outgoing particles are 
completely determined by four-momentum conservation: 

,/M 4 - 2M 2 {m{ + m 2 ) + (m 2 - m 2 ) 2 

~ P2 = 2M 

M 2 + m\ — m 2 , 
2M 

M 2 + mj- m 2 2 
2M 



E 1 = 
E2 = 



where M is the mass of the decaying meson and mi,m2 are the outgoing lepton masses. 
Several features of the kinematics are important. If one of the outgoing leptons has a mass 
close to the decaying meson, as the muon does in pion decay, then it takes most of the 
energy in the decay. This is the limit where mi — ► M above. In the pion decay chain 
then, the electron and three resulting neutrinos each carry about 1/4 of the pion energy. If, 
however, both outgoing leptons are light compared to the decaying meson, as is the case for 
kaon decay, the two outgoing leptons share equally in the energy. This is the limit where 
M ^> mi,TO2 above. The neutrinos which result from charged kaon decay, therefore, are 
higher in energy than those that result from pion decay. 
The energy spectrum in the lab frame is 

dN B 
dE~ v ~ 2pP lah 

where B is the branching ratio and Pi a b is the parent momentum in the lab frame. For the 
charged kaon decay, for example, this is 
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dN _ 0.634 

dE v ~ (1 - ml/M 2 K )P K 

To calculate the full flux of atmospheric neutrinos, one needs to solve coupled 
cascade equations for the nucleon and meson distributions that take into account interaction 
and decay lengths. The neutrino energy spectrum from each meson decay, like that for the 
kaon above, is then folded in with the meson distributions. This can be done with various 
analytical approximations (see chapters 3, 4, 6, and 7 of [12] for details) or it can be done 
with Monte Carlo techniques. The Monte Carlo techniques are most accurate and take into 
account the geomagnetic effect, which imposes an energy cutoff on primaries and bends the 
paths of charged secondaries, as well as the three-dimensional nature of the interactions. 
See [19, 20] for details on the canonical Bartol atmospheric neutrino flux which is used 
throughout this dissertation. The atmospheric muon and electron neutrino fluxes from the 
Bartol calculation are plotted in figure 1.2. Figure 1.3 illustrates their zenith dependence. 

1.3 Prompt Atmospheric Neutrinos 

At very high energies, charmed mesons like D and D° and charmed baryons 
like A+ will be produced in the primary cosmic ray interactions. These charmed particles 
have very short lifetimes and therefore very high critical energies (~ 10 7 GeV for the D's). 
Because the critical energies are so high, the energy spectrum of these so-called prompt 
neutrinos follows that of the primary cosmic ray spectrum and is isotropic with respect to 
zenith angle because the particles always decay before interacting. The decays also result 
in roughly equal numbers of muon and electron neutrinos. 

The only source of tau neutrinos from the atmosphere is a small component from 
the charmed, strange meson D~-. Calculations indicate that this prompt tau neutrino flux 
is exceedingly small, an order of magnitude lower than the prompt muon neutrino and 
electron neutrino flux [21]. 

The production of charmed mesons at cosmic ray energies involves hadronic physics 
which is not well-measured at accelerators. Flux models therefore vary over a large range. 
Figure 1.4 shows several prompt flux models superimposed on the conventional Bartol 
atmospheric neutrino fluxes. Figure 1.5 shows the zenith dependence of the prompt flux. 
Measuring the prompt atmospheric neutrino flux is a major goal of neutrino telescopes like 
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IceCube. From figure 1.4, it is clear that this should be easier to do with electron neutrinos, 
where the prompt component rises above the conventional component at almost an order 
of magnitude lower energy. 
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Figure 1.2: Atmospheric muon neutrino (left) and electron neutrino (right) fluxes from the 
Bartol model. The fluxes have been multiplied by E 2 and have units of GeV cm~ 2 s~ l sr _1 . 




Figure 1.3: Atmospheric muon neutrino (left) and electron neutrino (right) fluxes as a 
function of zenith and energy from the Bartol model. The fluxes have been multiplied by 
E 3 and have units of GeV 2 cm" 2 s" 1 sr" 1 . The fluxes are peaked at the horizon. 
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Figure 1.4: Conventional Bartol atmospheric neutrino flux (red) and prompt flux models 
(blue and green) for muon neutrinos (left) and electron neutrinos (right). The fluxes have 
been multiplied by E 2 and have units of GeV cm" 2 s" 1 sr -1 . Solid blue is the model from 
[22], green is the model from [21], and dashed blue is the model from [23]. 




Figure 1.5: Prompt flux for muon neutrinos (left) and electron neutrinos (right) as a function 
of zenith and energy for the model from [21]. The fluxes have been multiplied by E 3 and 
have units of GeV 2 cm" 2 s" 1 sr _1 . The fluxes are isotropic. 
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1.4 Neutrino Astrophysics 

As we have seen, cosmic rays are accelerated to enormously high energies in astro- 
physical sources. While evidence may be starting to confirm the long-held suspicion that 
galactic supernova remnants accelerate cosmic rays up to the knee (see section 1.5.3), the 
source of the highest energy cosmic rays is still unknown. However, the very existence of 
such high energy accelerators opens the door to the possibility of observing high energy 
astrophysical neutrinos. For this reason, neutrino astrophysicists have been designing and 
constructing large-volume neutrino telescopes in water and ice since the early 1990's. 

Astrophysical neutrino production proceeds along the same lines as the neutrino 
production in the atmosphere described above. The accelerated protons inside astrophysical 
accelerators should interact with ambient matter and photon fields to produce high energy 
astrophysical neutrinos via pion and kaon decay. This scenario is often described as a 
"cosmic beam dump" in analogy to the production of neutrino beams at particle accelerators 
on earth. 

The expected sources of these astrophysical neutrinos share several features. In 
general, they exhibit non-thermal emission, where the energy comes from the accretion of 
infalling matter or the gravitational collapse of a massive object. Sources which show very 
high energy (> TeV) gamma ray emission are also likely candidates. While gamma rays 
can be produced via purely leptonic mechanisms (accelerated electrons which synchrotron 
radiate and inverse Compton scatter), they can also come from hadronic interactions of 
accelerated protons via neutral pion decay: ir° — ► 77. Such a hadronic mechanism would 
also produce the charged pions that decay to neutrinos. 

A general argument can give us a sense of what these astrophysical accelerators 
might be. They can be classified by their size R and their magnetic field strength B. 
In order to be an efficient accelerator, an astrophysical source should be large enough to 
contain the orbits of the accelerated particles. This means that the source size should be 
larger than the gyro radius -Rgyro = -§ for a particle of energy E. The maximum energy 
achievable is then E = 'yBR where the gamma factor is included because the source may be 
boosted with respect to earth. Figure 1.6 shows the so-called Hillas plot of magnetic field 
strength as a function of source size needed to accelerate particles to a given energy. Also 
displayed are values for various source classes. 

If they can be detected from such sources, neutrinos would have several desirable 
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Figure 1.6: Hillas plot of magnetic field strength versus source size. The green line indicates 
what is necessary to accelerate an iron cosmic ray to 10 20 eV, and the red lines indicate 
what is needed to accelerate a proton to 10 20 -10 21 eV. From [24]. 
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features as astronomical messengers. They can probe very energetic, distant sources which 
are not detectable via other channels. High energy gamma rays from distant sources interact 
with background photons to produce electron-positron pairs and are hence attenuated. This 
occurs above a threshold given by 4£ , 7 i? 7 bg = (2m e ) 2 . For photons of energy 10-100 TeV 
interactions occur on the infrared background, and for PeV photons they occur on the cosmic 
microwave background. Because neutrinos interact only weakly, they are not attenuated or 
blocked by intervening photon fields or matter. In addition, neutrinos are not deflected by 
magnetic fields between the source and the earth and should therefore point directly back 
to their sources. High energy protons, on the other hand, are scrambled by magnetic fields 
and are also absorbed on the cosmic microwave background (the GZK effect). Because 
neutrinos have the ability to probe sources not readily observable in photons and protons, 
they may also reveal interesting surprises about the universe. 

This section will give a very brief introduction to the field of neutrino astrophysics. 
For good reviews of the subject, see [25, 26, 27, 28]. We'll begin with a tour of neutrino 
telescopes, continue with a discussion of Fermi acceleration and astrophysical neutrino pro- 
duction, and then review several likely sources of high energy astrophysical neutrinos. 

1.4.1 Prototype and Next-Generation Neutrino Telescopes 

To observe the low fluxes expected from astrophysical neutrino sources, neutrino 
telescopes need extremely large detector volumes, eventually reaching the cubic kilometer 
scale and beyond. The bulk of the existing and planned neutrino telescopes detect the 
Cherenkov light emitted by charged secondary particles produced in neutrino interactions. 
They therefore make use of large, naturally occurring volumes of transparent Cherenkov 
media — freshwater, sea water, and glacial ice. Figure 1.7 illustrates several of the completed 
detector arrays. 

The Baikal experiment [8, 30] is located 1.1 km below the surface of Lake Baikal 
in Siberia. It started operating in 1998 and was the first installed underwater neutrino 
telescope. Its NT200 configuration consisted of 8 "strings" , each of which had 24 optical 
sensors arranged in pairs to reduce backgrounds from radioactivity and bioluminescence. 
NT200+ added 3 additional strings around this dense inner core. In April 2008, Baikal 
deployed and successfully tested a prototype string employing new optical sensors and digital 
technology for full waveform readout [31]. The collaboration hopes that this technology will 
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Figure 1.7: ANTARES (left, copyright F. Montanet, CNRS/IN2P3 and UJF for Antares), 
AMANDA (center, from [29]), and Baikal (right, from [30]). 

form the basis for a future cubic kilometer scale detector. 

Other large-volume neutrino telescopes have been completed or are under con- 
struction in the Mediterranean. The ANTARES array [11, 32] lies 40 km off the coast 
of Toulon, France at a depth of 2500 m and was completed in May 2008. It consists of 
12 "lines", each with 75 optical sensors. They are arranged into "storeys", each of which 
has three sensors pointed downwards at an angle of 45°. Elsewhere in the Mediterranean, 
the NEMO and NESTOR collaborations are in the prototype and construction stages of 
building neutrino telescopes off the coasts of Sicily and Greece, respectively [33]. 

The ANTARES, NEMO, and NESTOR collaborations have joined together to 
pool their experience and resources to construct a cubic kilometer scale detector known as 
KM3NeT [34]. KM3NeT is still in the design, site selection, and planning stages. 

At the Amundsen-Scott South Pole Station in Antarctica, the AMANDA array 
took data from 1995 until 2009. AMANDA-B10 was composed of 10 strings with 302 
optical sensors embedded in the deep glacial ice. It was eventually upgraded to AMANDA- 
II, with a total of 19 strings and 677 optical sensors mostly deployed between 1500 m and 
2000 m below the surface of the ice. 

AMANDA served as the prototype array for IceCube, the first cubic kilometer 
neutrino telescope. At its completion, IceCube will have 86 strings, each with 60 optical 
modules capable of full waveform digitization in ice. IceCube is the primary subject of this 
dissertation and will be discussed in detail in chapter 3. 
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Because of the overwhelming background of downgoing muons from cosmic ray 
air showers, neutrino telescopes search for upgoing muons induced by neutrinos that have 
passed through the earth. For this reason, Baikal and the Mediterranean neutrino telescopes 
view the southern sky and have the galactic center in their fields of view. AMANDA and 
IceCube, on the other hand, view the northern sky. From the point of view of coverage, it 
is therefore desirable to have cubic kilometer scale neutrino telescopes in both hemispheres. 

1.5 Neutrino Production in Astrophysical Sources 

Astrophysical neutrinos are believed to come from the decay of mesons produced 
in the interactions of accelerated protons with ambient matter and radiation fields in and 
around astrophysical accelerators. The accelerated protons can be linked with the observed 
highest energy cosmic rays, or they can occur in hidden sources where the protons and 
photons don't make it out of the source but the neutrinos do. The mechanism of Fermi 
acceleration is thought to be responsible for accelerating these protons and for giving them 
their characteristic power law spectrum. 

1.5.1 Fermi Acceleration 

The essential idea of Fermi acceleration is that a moving plasma can transfer 
bulk kinetic energy to individual charged particles. The charged particles diffuse through 
turbulent magnetic fields in the moving plasma, scattering elastically in encounters with 
the magnetic field. 

To get the basic contours of Fermi acceleration, we follow the simple argument 
from chapter 11 of [12]. At each encounter in the magnetized plasma, the charged particle 
gains an amount of energy proportional to its energy: 

AE = £E 

After n encounters, the energy is 

E n = E Q (l + O n 

In order to reach an energy E, then, the particle must undergo a number of encounters 
given by 
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71 = — 

ln(l+0 

Next, we assume that the particle has some probability P esc of escaping the accelerator at 
each encounter. The probability of a particle remaining after n encounters is then (1— P esc ) n . 
In order to reach an energy E or greater, the particle has to survive n or more encounters 
where n is given by the equation above. So we can write 

oo 

N{> E) CC ^(1-Pesc) 
m=n 

Taking the natural log of both sides we have 



(1 - P a 



In N{> E) oc In ( — J + n ln(l - P csc ) 



PcscJ \E J ln(l+{) 



Exponentiating again we have 



N(> E) oc J- 



where 

hi 



1— Pesc I Pcsc _ 1 ^cycle 



7 ln(l + ~ e e X ^csc 
and T cyc i e and T esc are the characteristic acceleration-cycle and escape times. Fermi accel- 
eration naturally leads to the power law spectrum that we were trying to reproduce for the 
observed cosmic rays. 

To determine the actual power law index, the fractional energy gain at each en- 
counter must be calculated by considering the scattering of particles from different magnetic 
field configurations. Fermi acceleration can be divided into two types. In first order Fermi 
acceleration, particles are accelerated by a large, planar shock front which is moving at a 
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velocity v. The fractional energy gain is proportional to (3 = v/c. In second order Fermi 
acceleration, particles are accelerated by a moving gas cloud, and the fractional energy gain 
is proportional to (3 2 . First order Fermi acceleration at shock fronts in supernova blast 
waves is thought to be responsible for accelerating the bulk of the cosmic rays up to the 
knee, and it is also a possible mechanism for accelerating very high energy protons and 
electrons in gamma ray bursts and active galactic nuclei, which will be discussed below. 

1.5.2 Proton-Proton and Proton-Photon Interactions 

To produce neutrinos, these Fermi accelerated protons interact with ambient mat- 
ter or radiation fields in and around the source. These interactions can be divided into two 
types: proton-proton (or simply p-p interactions) and proton-photon (p-7) interactions. 
The products of these reactions are 



If the density of the astrophysical source is sufficiently low, these mesons decay in the 
familiar chains described above before they have a chance to interact. 

Consider p-p interactions first. The "volume emissivity" of pions in units of 
GeV -1 s" 1 cm -3 can be expressed as [26]: 



where n t is the density of protons in the target, N n is the produced pion multiplicity, 
dcr p p(E n , E) I dE is the differential inclusive cross section for producing N n pions of energy 
E w from a proton of energy E, and dN/dE is the proton spectrum in units of GeV -1 cm" 3 . 
We can simply estimate the cross section as da pp (E n , E)/dE ~ a pp SlE^ — K p E/N n ] where 
a pp ~ 30 mb is the proton-proton inelastic cross section at TeV lab energies, k p ~ 0.4 is the 
inelasticity, and N n ~ 15. This says that the pions share the transferred energy equally. 
Then the integral above becomes 



The neutrino volume emissivity follows from this and from kinematics. Roughly 1 /3 of the 
proton energy goes into each of the charged and neutral pions. In a thin source, the pions 



P + P>P + 7 — ► 7T 



,ir ± ,K ± + X 
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decay before interacting and three neutrinos are produced, each with roughly 1/4 of the 
parent pion energy (in dense sources pion and kaon energy loss and interaction must be 
taken into account). So we can write 

Ql p ~ 3 x 1 x 2 -Q™{AE v ) 
The neutrino flux then follows by integrating over the spatial extent of the source: 




where & u has units of GeV s cm 

For p-7 interactions, a similar calculation can be carried out for a target density 
of photons with a blackbody spectrum or a power law spectrum (which might come, for 
example, from synchrotron radiation). We need to impose the threshold condition for pion 
production in the intgral and to use a pl « 120 fih well above threshold, k p « 0.3, and 

In p-7 interactions, resonant production of the A is particularly important. This 

9 9 

occurs above an energy threshold given by E p E~ f = — ^ — -. Flux calculations generally use 
the fact that < x p ^^ >« i is the average fractional energy transfer from the proton to the 
pion produced through the delta resonance. 

1.5.3 Galactic Sources 

Several galactic candidates for neutrino emission have been considered in the liter- 
ature. Pulsars are rapidly spinning neutron stars with strong magnetic fields. Electrostatic 
acceleration of particles from the polar surface regions of neutron stars has been discussed 
as a source of high energy cosmic rays. If such acceleration occurs, the protons can interact 
with x-rays, also produced by the pulsar, to produce pions and neutrinos [35]. 

Microquasars are several solar mass black holes accreting matter from a binary 
companion. They show hard x-ray emission as well as relativistic jets. Models of neutrino 
emission include p-7 interactions in the jets [36] and p-p interactions with matter in the 
accretion disk or in the companion star [37]. 

Perhaps the most exciting galactic candidates are the shell type supernova rem- 
nants, which are thought to be the main acceleration sites of galactic cosmic rays. Recent 
gamma ray observations by the air Cherenkov telescope H.E.S.S. of the shell type supernova 
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remnants RX J1713. 7-3946 and RX J0852. 0-4622 show strong TeV gamma ray emission. 
Modeling of the multiwavelength spectrum from the radio through x-rays and gamma rays 
suggests that the gamma ray emission is dominated by the decay of ir° produced in the 
interactions of shock-accelerated protons [38, 39, 40]. Observations of neutrinos would be a 
smoking gun for the acceleration of high energy cosmic rays in these shell type supernova 
remnants [41]. 

1.5.4 Extragalactic Sources 
Active Galactic Nuclei 

Active galactic nuclei (AGN) are million to billion solar mass black holes situated 
at the centers of galaxies. They are the brightest objects in the universe, with some radiating 
as much power as all the stars in the Milky Way from a region smaller than the solar system 
[26]. The energy source for these objects is the gravitational energy released from accreting 
matter. Many AGN exhibit relativistic jets, TeV gamma ray emission, and short timescale 
flares. When the AGN jet points towards the earth, it's known as a blazar. 

Shocks in these jets are possible sites for the acceleration of protons. The protons 
can interact with matter in the jet itself or with the thermal ultraviolet photons from the 
heated infalling matter to produce neutrinos. 

Gamma Ray Bursts 

Gamma ray bursts are extremely bright, cataclysmic explosions distributed over 
cosmological distances that last on the order of % ~ 10~ 3 -10 3 seconds. The non-thermal 
gamma ray emission is characterized by broken power laws and can extend into the GeV 
energy range. Long gamma ray bursts (tj, > 2 s) are associated with the collapse of very 
massive stars to black holes. Short gamma ray bursts (tb < 2 s) are thought to arise from 
the mergers of neutron star binaries or neutron star-black hole binaries. In both cases, 
the newly formed black holes can accrete matter and form highly beamed jets with bulk 
Lorentz factors of ~ 100. See [42] for a good review of the observational and theoretical 
issues surrounding GRB's. 

The jet is an optically thick e^,"/ "fireball". The gamma ray photons can result 
from synchrotron radiation and inverse Compton scattering from electrons which are Fermi 
accelerated at internal shocks in the jet. In addition, the jet may contain accelerated baryons 
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(i.e. protons). Accelerated protons can interact with these MeV synchrotron photons to 
photoproduce mesons and neutrinos through the delta resonance described above. This 
leads to neutrino energies in the PeV range. Also, if the jet had to burrow through the 
collapsing precursor star on its way out, proton-proton and proton-photon interactions can 
lead to precursor neutrinos 10-100 s before the observed gamma rays. 

Choked Supernovae 

Gamma ray bursts have been associated with type Ib/c supernovae. While the 
jets of these gamma ray bursts have Lorentz factors ~ 100, it has been suggested by [43] 
that mildly relativistic jets with Lorentz factors on the order of a few may be a more general 
feature of supernova collapse. These mildly relativistic jets may not make it all the way 
through the stellar envelope, as in a GRB. However, such "choked bursts" will still produce 
neutrinos without an observable gamma ray signature [44]. 

1.6 Oscillations and the Astrophysical Flavor Ratio 

We now have abundant evidence for the phenomenon of neutrino oscillations from 
atmospheric, solar, and accelerator neutrinos. That is to say, a neutrino which is produced 
in a weak interaction in a state of definite flavor has a nonzero probability of being detected 
as a neutrino of a different flavor in a detector which is situated some distance from the 
source. This is only possible if neutrinos have mass. Oscillations will alter the flavor ratios 
of astrophysical neutrinos detected on earth. 

We can think of two sets of eigenstates for neutrinos — the flavor basis and the 
mass basis. These sets of eigenstates are rotated with respect to one other by a unitary 
mixing matrix U* [17]: 

K >=^2km > 

i 

where the index a refers to the flavor eigenstates and the index i refers to the mass eigen- 
states. The mixing matrix U is parameterized by three angles and can be written as 
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where Sjj = sin Oij and Qj = cos 9y, the 0« are the mixing angles, 5 is a possible CP- violating 
phase, and a±, a 2 are possible Majorana phases. 

Each mass eigenstate picks up a phase exp(— irriiTi) during propagation, where mi 
is the mass of the mass eigenstate and Tj is proper time. The probability for v a to oscillate 
into V13 can be found by multiplying each term in the sum for \v a > with its phase factor and 
then projecting the state < vp\ onto it. For details see [17]. The formula for the oscillation 
probability is 

P(u a ->■ vp) =5 a/3 

- ^^{UyjpiU^) sin 2 [1.27Am2.(L/ J E)] 

i>j 
i>j 



where Amfj = m 2 — mj is in eV, L is in km, and E is in GeV. 

For astrophysical neutrinos, oscillations mean that the flavor ratio at a detector 
on earth will be different from the flavor ratio generated at the astrophysical source. From 
the meson decay chains described above, we expect a generic flavor ratio at the source of 

(z^ : v e : z^r) source « 2 : 1 : 

The effect of oscillations over the extremely long baselines of astrophysical sources is to 
transform this ratio to 
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{v„ : v e : ^ r ) earth 

That is, we expect an equal ratio of fluxes at our detector on earth. For details of the 
calculation, see Appendix I. 

This has an important implication for astrophysical neutrino detection. Even 
though we don't expect astrophysical sources to produce tau neutrinos in appreciable num- 
bers, they should be present in an astrophysical neutrino beam because of oscillations. Since 
the atmospheric tau neutrino flux is small (no conventional contribution but an exceedingly 
small contribution from prompt charmed mesons), the detection of a tau neutrino would be 
an almost smoking gun for an astrophysical source. 

In addition, by detecting the different interaction signatures from different neu- 
trino flavors from an astrophysical neutrino source, neutrino telescopes can hope to probe 
oscillations over extremely long baselines for effects of exotic physics. This is described in 
more detail in the next chapter. 
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Chapter 2 

Principles of Neutrino Detection 

Neutrinos are often described as "ghost particles" because they are so difficult 
to detect. In fact, several of the big names of twentieth-century physics thought neutrino 
detection would never happen. "It is therefore absolutely impossible to observe processes of 
this kind with the neutrinos created in nuclear transformations," Hans Bethe and Rudolph 
Peierls wrote of the possibility of observing neutrinos via inverse beta decay [45]. Never- 
theless, over the course of the last half-century, physicists and astrophysicists have become 
increasingly adept at building the large-scale detectors needed to observe inverse beta decay 
and other types of neutrino interactions with reasonably high statistics. 

This chapter will discuss the underlying physics that goes into such detection. 
First, we'll give a brief description of Cherenkov light, the basic physical effect that un- 
derlies many different neutrino experiments, including IceCube. Next, we'll discuss the 
fundamental interactions of neutrinos and how they manifest themselves at different energy 
scales. To understand the rates of each interaction we'll discuss the neutrino cross sections. 
Finally, we'll go on to talk about charged particle energy loss and conclude with one of the 
particular signatures of very high energy neutrinos — electromagnetic and hadronic particle 
showers or "cascades". Cascade detection is the primary subject of this dissertation. 

2.1 Cherenkov Light 

Cherenkov light is produced when a charged particle travels through a dielectric 
material faster than light can pass through the same material. It is the optical analogue of 
a sonic boom and is responsible for the faint blue glow that emanates from a nuclear reactor 
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core. The wavefront of Cherenkov light forms a cone whose axis is the particle's trajectory. 
The wave vectors make an angle 9 C with the particle's track: 

cos C = — 
pn 

for > — , where j3 = v/c and n is the index of refraction of the material. This last 
inequality expresses the condition that the particle's velocity exceed the speed of light 
in the material and leads to a momentum threshold below which a particle cannot emit 
Cherenkov radiation: 



In water or ice with n = 1.33, the momentum threshold is 0.58 MeV/c for an electron, 120 
MeV/c for a muon, 159 MeV/c for a pion, 563 MeV/c for a kaon, and 1070 MeV/c for a 
proton. For (3 ~ 1, the Cherenkov angle in water and ice is 9 C ~ 41°. 

The amount of Cherenkov light emitted along a particle's trajectory is given by 
the Frank- Tamm formula: 



where a is the fine structure constant and z is the particle's charge in units of the electron 
charge. In ice, an electron or muon traveling with (3 ~ 1 emits around 330 photons per 
centimeter in the 300 to 600 nm range of interest for IceCube's optical sensors. 

2.2 Neutrino Interactions and Cross Sections 

Depending on the source, neutrino energies vary over some nine orders of mag- 
nitude, from MeV-scale reactor, solar, and supernova neutrinos to GeV-scale atmospheric 
and accelerator neutrinos into the TeV- and PeV-scale astrophysical neutrinos that IceCube 
hopes to detect. For different neutrino energies and targets, various classes of reaction are 
possible between the neutrino and the detector medium. 





At their most basic level, neutrino-quark reactions are of two fundamental types: 
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v\ + q — > I + q' Charged-Current (CC) 
vi + q — ¥ v\ + q Neutral-Current (NC) 

where I stands for the lepton flavor and q, q' are quark species. The charged-current 
interaction is mediated via exchange of the charged W boson, and the neutral-current 
interaction is mediated via exchange of the neutral Z boson. These fundamental reactions 
manifest themselves in different ways depending on the neutrino energy. For a pedagogical 
introduction, see [46]. 

2.2.1 The Cross Section 

In order to predict how many events of each class we expect in a detector, we need 
to know the cross sections for neutrino reactions. The number of events will be given by 

N V {E) ~ <b v {E) x a v {E) x iV targct x T 

where <& V {E) is the flux of neutrinos of energy E, a u (E) is the cross section, N ta _ Tget is the 
number of targets (free protons, deuterons, bound nuclei like carbon and oxygen, electrons, 
etc.), and T is the observation time. 

The next few sections will give a brief introduction to the neutrino reactions that 
are relevant at different energy scales as well as their cross sections. Particular attention 
will be paid to deep inelastic scattering, the regime most important for IceCube. 

2.2.2 Low Energy 

At the lowest energies relevant for reactor, solar, and supernova neutrinos (< 100 
MeV), the important processes are elastic and quasi-elastic scattering. Since there are no 
free neutrons in materials, the basic free nucleon reaction is inverse beta decay: 

v e + P — ► e+ + n 

Inverse beta decay has a threshold around 1.8 MeV and a distinct signature from the 
delayed-coincidence of the annihilation of the positron and the neutron capture. It was 
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the basis of the first neutrino detection experiments of Reines and Cowan [5, 6] and the 
reactor neutrino experiment KamLAND [2]. It is also the basis of upcoming reactor neutrino 
oscillation experiments like Double Chooz [47] and Daya Bay [48] that are trying to measure 
#13, the last unknown oscillation parameter. 

While inverse beta decay dominates the cross section at these energies (a ~ 
10 -41 cm 2 for 10 MeV neutrinos on water), other reactions are possible. Neutrinos can 
elastically scatter off electrons in the detector medium. The cross sections are smaller 
(a ~ 10™ 43 — 10 -44 cm 2 for 10 MeV neutrinos on water), but the recoiling electron pre- 
serves directional information from the incoming neutrino. This reaction was used to make 
the famous SuperK images of the sun in neutrinos [49] . 

Neutrinos can also interact with bound nucleons, though the energy thresholds 
are higher and the cross sections lie between the inverse beta decay cross section and the 
electron scattering cross section. For the deuterons in heavy water in experiments like SNO 
[3] the following reactions are important: 



v e + d — > e +p + p 
v e + d — > e + + n + n 



For large water detectors reactions on oxygen need to be taken into account: 



u e + 16 O -» e~ + 16 F 
u e + 16 O -» e+ + 16 N 



And in liquid scintillator there are also contributions from carbon: 



v e + 12 C -> e~ + 12 N 



D e + 12 C e+ + 12 B 



Radiochemical neutrino reactions on chlorine provided the first detection of solar v e by Ray 
Davis and colleagues [50]. 
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2.2.3 Intermediate Energy 

At energies on the order of a GeV, quasi-elastic scattering is still important, and 
there is now enough energy available to produce an outgoing muon in z/„ interactions as 
well: 

z/ M + n ->■ yT + p 
v e + n — > e~ + p 

plus the corresponding reactions for v e and Pa- 
li the neutrino has enough energy, it can also excite resonances (the most important 
of which is the A(1232) resonance) which lead to the production of charged and neutral 
pions: 

N* -> 7T + N' 

Production of pions can also proceed via coherent interactions with entire nuclei. Pion 
production is particularly important for accelerator experiments like NOz^A [51] and T2K 
[52] which are searching for the appearance of v e in a beam of v^. Final state 7r° produced 
in Vn interactions decay to photon pairs which can mimic the electrons produced in v e 
interactions. 

As energy increases in this intermediate regime, the total cross sections are linear 
with neutrino energy and are roughly [17]: 



a (uN) = 0.677 x 10~ 38 cm 2 x E„/(GeV) 
a[pN) = 0.334 x 10" 38 cm 2 x E„/(GeV) 
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2.2.4 High and Ultra-High Energy 

At 100's of GeV and above, we enter the deep-inelastic scattering (DIS) regime, 
where the neutrino has enough energy to scatter from quarks themselves, breaking up the 
nucleon in a hadronic shower. To write down the cross section for the charged-current DIS 
process i>i + N — > l + X, we introduce the conventional kinematic variables. The four-vectors 
are given by the following: 

incident neutrino: p = {E v ,p v ) nucleon: px = (Mjv,0) 

outgoing lepton: p' = (Ei,pi) outgoing hadronic final state: px = (Ex,Px) 



First, we define the total center-of-mass energy: 



s = (p + Pn) 2 = 2M N E U + M 2 N ^ 2M N E V 

Next, we define q as the four-momentum transfer to the W boson: q = p — p' = {y,q). We 
define Q 2 as the negative magnitude of the four-momentum transfer: 



Q 2 = -q 2 = -(p- p') 2 = 2{E V E X - p v ■ pi) -ml-m 2 ^ ±E V E X sin 2 -j 

where the last approximation holds when we can ignore the lepton masses compared to the 
energies and Oi is the angle of the outgoing lepton to the incoming neutrino. Next, we define 
the energy transfer to the W boson in the nucleon rest frame: v = = E v — E[. Finally, 
we can define the Bjorken scaling variables x and y: 



-q 2 _ Q 2 
1q ■ px 2Mxv 

q-PN _ ^ Q 2 Q 2 

p-p N E v 2M N E u x (s-M%)x 



The Bjorken x variable varies from zero to one. It's zero when Ei — > 0, that is all of the 
neutrino four-momentum is transferred to the W boson. It's one when E\ — > E v , that is 
all of the neutrino four-momentum goes into the outgoing lepton. x is also the fraction 
of nucleon momentum carried by the struck quark. This can be seen via the following 
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argument. Let £ be the fraction of nucleon momentum carried by the struck quark. At 
the vertex where the W boson interacts with the struck quark and produces an outgoing 
hadron we have 



= p\ = (&N + q? = q 2 + 2& ■ PN = ~Q 2 + 2£q ■ p N £ = = x 

2q ■ p N 

The Bjorken y variable is also called the inelasticity and is the relative energy transfer from 
the incoming neutrino through the W boson and into the nucleon in its rest frame. 

We're now in a position to write down the charged-current DIS cross section on 
an isoscalar target [53]: 

[x q (x,Q 2 ) + xq(x,Q 2 )(l-yf] 



dxdy vr (1 + Q 2 /M 2 W ) 2 

where the Fermi constant Gf = 1.16632 x 10~ 5 GeV~ 2 , M is the nucleon mass, and the 
quark distribution functions q(x,Q 2 ) and q(x,Q 2 ) are linear combinations of the valence 
and sea quark distributions for the different quark types. To obtain the neutral-current DIS 
cross section, we replace Mw with Mz and take another linear combination of the valence 
and sea quark distributions for the quark types. 

The cross sections are shown in figure 2.1. They rise linearly with energy up to 
several TeV and then decrease in slope because of the Q 2 /M 2 \y term from the propagator. 
Extrapolations must be made to carry the quark distributions into unmeasured regions of 
the small x, large Q 2 phase space. IceCube uses the CTEQ5 (Coordinated Theoretical- 
Experimental Project on QCD) distribution functions [54, 55] for calculating cross sections. 

One more process must be taken into account in the high energy neutrino cross 
section. At E 9e « M 2 w /2m e = 6.3 PeV, v e can resonantly scatter off electrons to produce 
the W~. This is the so-called Glashow resonance [56]. The cross section for this resonance 
is several hundred times larger than the charged-current cross section at this energy. The 
Glashow resonance is important for high energy astrophysical neutrino searches but doesn't 
come into play at the energies of atmospheric neutrinos. 

2.3 Charged Particle Energy Loss in Matter 



The interactions described above are the fundamental processes that allow for 
neutrino detection. In an actual experiment, one makes observations of the products of 
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log(Energy) [GeV] 

Figure 2.1: Neutrino-nucleon cross sections using the CTEQ5 distribution functions. 

each reaction — in particular, the muons, electrons, and recoiling nuclei. 

Charged particles traversing matter lose energy via ionization and radiation pro- 
cesses. For heavy particles other than electrons, the ionization loss is given by the Bethe- 
Bloch equation [57]: 



dE\ UitN A z 2 e A \ (Z_ 
dx )\on V m eV 2 



( 2m e v 2 j 2 



In 



where m e is the electron mass, v and z are the velocity and charge of the incident particle, 
Na is Avogadro's number, and / ~ WZ eV is the material's mean ionization potential. This 
expression goes as ~ v ~ 2 at low velocities, reaches a minimum at around 3Mc 2 where M is 
the charged particle mass, and then rises logarithmically with energy, eventually flattening 
out at a value around 2 MeV g _1 cm 2 . 

At higher energies, radiative losses like bremsstrahlung, pair production, and pho- 
tonuclear interactions start to kick in. The energy at which the ionization losses equal these 
radiative losses is the critical energy, which for a muon in a solid or liquid is given by [17]: 



E, 



fIC 



5700 
(Z + 1.47) 



0.838 



GeV 
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For water and ice this is around 1 TeV. The main background for all IceCube analyses 
consists of high energy muons and muon bundles produced by cosmic rays interacting in 
the atmosphere. These air shower muons are generally above the critical energy and so lose 
most of their energy radiatively. Most of a muon's Cherenkov light will therefore come from 
the secondary particles produced in these radiative losses. 

For high energy electrons, energy loss is dominated by bremsstrahlung. This energy 
loss can be expressed as 

dE\ E 

where Xq is the radiation length, the distance over which an electron loses on average all 
but a fraction 1/e of its energy. An approximate formula for the radiation length of an 
electron in units of g cm -2 is given by [17]: 

7WAA 



Z(Z + l)1n(287/VZ) 

which is around 37 g cm -2 or 41 cm for an electron in ice (density 0.9 g cm" 3 at 0° C). 

The energy at which the ionization losses equal the radiation losses is known as 
the critical energy E c and, for an electron in a solid or liquid, is given by [17]: 

E c = - 610 MeV 
(Z + 1.24) 

This is around 72 MeV for an electron in ice. 

The radiation length also controls the distance scale over which a photon converts 
to an electron-positron pair. On average, a photon pair-produces over a distance A pa i r = 
(9/7)X . 



2.4 Electromagnetic Showers or "Cascades" 

When a high energy electron passes through a material, it will emit a high energy 
photon through bremsstrahlung. This photon can convert to an electron-positron pair 
which, in turn, can emit photons of their own. This process goes on and on, forming 
an electromagnetic particle shower or "cascade" of high energy photons, electrons, and 
positrons. The high energy electrons that are produced in v e DIS interactions in a detector 
like IceCube manifest themselves as electromagnetic cascades. 
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To understand the main scaling features of an electromagnetic cascade, we can 
consider an extremely simple branching model [57]. We start with an electron of energy 
Eq which, after one radiation length Xq, radiates one bremsstrahlung photon of energy 
Eq/2 (in reality, this fraction varies, but in our simple model we'll always take it to be 
1/2). In the next radiation length, this photon converts to an electron-positron pair, each 
of energy Eq/A. The original electron also radiates another photon of energy Eq/A in this 
radiation length. This process continues to subdivide the energy until the particle energies 
fall below the critical energy E c , at which point all of the particles simply dissipate their 
energy through ionization. 

After t radiation lengths, each particle has energy E(t) = Eq/2 1 . We can define 
the "shower maximum" as the point at which the shower contains the most particles. In 
our simple model, this will occur at the point when the particles reach the critical energy 
at a depth i max given by 

E c = £ /2* max => t m ax = ln(£ /£e)/hi2 

So we see that the depth of shower maximum scales only logarithmically with the energy 
of the incident particle. The number of particles at the shower maximum is given by 

iVmax = 2* max = e' max 1112 = E /E c 

The number of particles at shower maximum, however, scales linearly with energy. The 
total track length of charged particles is given by 

where the factor of 2/3 roughly accounts for the charged particle fraction in our simple 
model. The total track length is again proportional to the incident particle energy. The 
total track length of a cascade is an important quantity because it will tell us the total 
amount of Cherenkov light emitted by the cascade. 

Putting in numerical values, we obtain a total track length of 13.4 Xo/GeV or 
5.5 m/GeV. Realistic cascade simulations have determined the actual total track length to 
be 5.9 m/GeV [58]. However, since the amount of Cherenkov light emitted by a particle 
depends on its velocity, we define an effective total track length to take into account the 
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velocity distribution of particles within the cascade. From simulations, this effective total 
track length is 5.2 m/GeV [58]. 

The longitudinal profile of an electromagnetic cascade can be fit by the following 
functional form [17]: 

dE Jtb^e-* 

H ~ E ° b r(o) 

where t is the depth in units of radiation lengths, Eq is the electron energy, and a and b are 
constants which, for ice, were determined by fits to a GEANT simulation [58, 59]: 



a = 2.03 + 0.604 • In 



Eq 
GeV 



b = 0.633 



The shower maximum is given by t max = (o — !)/&■ Almost all of the energy deposition 
occurs within 20 radiation lengths. 

In realistic simulations of an electromagnetic cascade, Coulomb scattering of the 
electrons and positrons as they pass through the medium must be taken into account. This 
gives a lateral spread to a shower, which is governed by the Moliere radius Rm [17]: 

Rm = ( 21 ' 2 ^ MeV ) = 0.0265(Z + 1.2)X 

In ice, we have Rm ~ 8.3 g cm -2 or 9.2 cm. 90% of the shower energy is deposited within 
a cylinder of radius Rm-, and 99% is deposited within a cylinder of radius 3.5 Rm- 

To summarize then, an electromagnetic cascade in IceCube deposits the bulk of 
its energy in a thin cylinder ~ 30 cm in radius and ~ 5 m in length. 

2.5 Hadronic Showers or "Cascades" 

Energetic hadrons moving through a material result in a shower of baryons and 
mesons. These are known as hadronic showers or cascades. The recoiling nuclei in high 
energy neutrino DIS will produce a jet of fragmented particles, each of which will shower. 

In general, a hadronic cascade produces less Cherenkov light than an electromag- 
netic cascade of the same energy. This happens for several reasons. First, some of the 
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hadronic cascade energy will go into neutrons, which produce no Cherenkov light. Second, 
some of the energy is dissipated in the large nuclear binding energies involved. Finally, 
pions, kaons, and protons have higher Cherenkov thresholds than electrons and positrons 
do [58]. 

However, the relative light yield between hadronic and electromagnetic cascades is 
a function of energy. The hadronic cascade will produce 7r°'s, which decay to photon pairs. 
These photons produce an electromagnetic cascade and represent a "one-way street" [60] 
by which energy leaves the hadronic sector for the electromagnetic sector. As the hadron 
energy increases, more 7r°'s are produced and the cascade becomes more electromagnetic. 

Defining the relative total track lengths for a hadronic and electromagnetic cascade 

as [61] 

7-, -^hadronic 

= T 

we can write 

F = Fcm + (1 — Fcm) • fo 

where F ern is the electromagnetic part of the cascade fed by 7r°'s, 1 — F cm is the purely 
hadronic part of the cascade, and /o represents the relative total track length for a pure 
hadronic cascade. Following [60], we model the energy dependence of the electromagnetic 
fraction F cm as a power law 

F cm = 1 - (E/E )- m 

where E is the hadron energy and Eq and m are fit parameters determined individually 
from Monte Carlo simulations of injected p's, n's, 7r's, and K J s [61]. 

In a realistic deep inelastic scattering event, the recoiling nucleus will fragment 
into a jet of particles, each of which will individually shower. To determine Fx, the total 
track length fraction for this recoiling nucleus at the neutrino-interaction vertex, simula- 
tions were performed by the AMANDA collaboration to determine the composition of the 
fragmented particles. Fx was obtained by superimposing the total track length ratios from 
the individual components of the jet (p's, n's, 7r's, iT's, etc.) [61]. The mean and RMS 
values of Fx as a function of the hadronic cascade energy are illustrated in figure 2.2. 
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Figure 2.2: Mean and RMS of the total track length ratio F for hadronic cascades. 

2.6 Event Signatures in IceCube 

We are now in a position to put everything together to discuss the main neutrino 
event signatures in the IceCube detector, which detects Cherenkov light from the charged 
secondary particles from neutrino-nucleon interactions. The detector itself will be described 
in detail in the next chapter. 

Figure 2.3 shows IceCube event- viewer snapshots of the three main neutrino event 
signatures. The IceCube event-viewer depicts the amount of detected Cherenkov light in a 
sensor. The size of a balloon indicates the amount of detected light, and the colors indicate 
timing — from the earliest hits in red to the latest hits in blue. 

The left panel is a muon track in IceCube. This particular muon is from a cos- 
mic ray air shower background event and is moving downward through the detector, but a 
neutrino-induced muon would have the same appearance moving upward through the de- 
tector. Cherenkov light is emitted from the track itself as well as from the radiative energy 
losses along the track. 

IceCube observes muon neutrinos by looking for muons passing upwards through 
the detector. An upward going muon track must have come from a neutrino, since only a 
neutrino could have traversed the interior of the earth to interact in the rock or ice below 
the detector. Cosmic ray air shower muons on the other side of the planet are all absorbed. 

The middle panel shows a simulation of an electron neutrino cascade event. The 
Cherenkov light comes from the relatively local energy deposition of the electromagnetic 
cascade from the outgoing electron plus the hadronic cascade at the interaction vertex. 
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Figure 2.3: Different neutrino flavors have different signatures in IceCube. A long, straight 
muon track is the hallmark of a muon neutrino (left). A point like deposition of light from 
a cascade (center) is produced by electron neutrinos and neutral-current interactions of any 
neutrino flavor. Tau neutrinos have a distinctive double bang signature (right). 

These two separate cascades cannot be resolved individually. In a sparsely-instrumented 
detector like IceCube, cascades look effectively like point sources of Cherenkov light. 

Several other classes of neutrino events also have the cascade topology. First, 
the hadronic cascades from neutral-current interactions of all neutrino flavors appear as 
cascades. Second, a muon neutrino which has a charged-current interaction inside of the 
detector will emit light from the hadronic cascade at its interaction vertex. If the inelasticity 
is large, the hadronic cascade can dominate the light output from the relatively dim outgoing 
muon. These are called "starting" muon neutrino events. Finally, low energy tau neutrino 
charged-current interactions appear as cascades. 

The right panel shows a simulated tau neutrino event. Tau neutrino charged- 
current interactions produce a hadronic cascade at the interaction vertex plus an outgoing 
tau lepton. The tau lifetime is 2.9 x 10~ 13 s, but if it has enough energy it can travel an 
appreciable distance before decaying (a 1 PeV tau will travel around 50 m before decaying). 
Most of the time the tau decay results in another cascade. This topology, where the two 
cascades are separated by the long tau track, is known as a "double bang". As the tau 
neutrino energy decreases, the two cascades move closer together and eventually become 
indistinguishable from a single cascade. 

As we discussed in section 1.6, tau neutrinos have one distinct advantage over 
muon and electron neutrinos: they are not produced in the atmosphere (except for an 
exceedingly small contribution from prompt charmed meson decay [21]). A detection of tau 
neutrinos would be a smoking gun for an extraterrestrial neutrino source. 
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2.7 Why Search for Cascades? 

There are several reasons why we would like to be able to identify the cascade 
event signature in a high energy neutrino telescope like IceCube. First and foremost, in 
the search for an extremely small signal, one wants to be as efficient as possible. Neutrino 
oscillations should cause any astrophysical neutrino flux to reach earth in a flavor ratio of 
~ 1 : 1 : 1 (see section 1.6), and each neutrino flavor can interact with the cascade signature. 
All electron neutrinos appear as cascades, muon and tau neutrino neutral-current events 
manifest themselves as hadronic cascades, and low energy tau neutrino charged-current 
events look like single cascades. So an appreciable fraction of any astrophysical signal will 
make itself known in the cascade channel. Searching for cascades in conjunction with muon 
tracks should improve our sensitivity to any small astrophysical signal. 

If astrophysical neutrinos are detected, a natural and important next task is to 
measure their flavor ratio. In order to tease out the flavor ratio, one needs to be able to 
measure muon tracks, cascades, and tau signatures. While this is challenging because of 
the challenges of cascade and tau detection, the payoff is large. Several factors can alter 
the flavor ratio from the expected 1:1:1 ratio at the detector. New physics like neutrino 
decay, CPT violation, sterile neutrinos, and pseudo-Dirac neutrinos are all expected to alter 
the astrophysical neutrino flavor ratio [62, 63, 64, 65, 66, 67]. In addition, the flavor ratio 
is sensitive to the physics and the environment of the astrophysical neutrino source. The 
charged pions, kaons, and muons that decay to neutrinos can lose energy in the source 
due to interactions with radiation and magnetic fields. This alters the flavor ratio [68]. 
The production and subsequent decay of short-lived charmed mesons like D ± and D° m 
astrophysical sources will also alter the flavor ratio [69]. So by measuring cascades and 
the flavor ratio of astrophysical sources, we can probe new physics and the physics of 
astrophysical sources. 

Next, cascades are useful for astrophysical searches over the full sky. In this way, 
cascade searches are quite complementary to muon searches. Muon neutrino searches have 
very good angular resolution. The muon track follows the incoming neutrino direction 
with a deviation given by 0.7° /(-E^/TeV) ' 7 , and the track directional reconstruction has 
a resolution ~ 1°. In practice though, the enormous background of cosmic ray air shower 
muons makes it extremely difficult to search for astrophysical muon neutrinos from the 
southern sky. Cascades, on the other hand, have poor angular resolution (see chapter 
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4). However, if the cascade event topology can be separated from the cosmic ray muon 
background, then it can be done over the full sky. Already, a search for neutrino-induced 
cascades from gamma ray bursts over the full sky has been conducted with AMANDA [70]. 

Finally, cascades are particularly powerful in searches that look for a diffuse flux. 
This diffuse flux could come from the undetected prompt component of the atmospheric 
neutrino flux (see section 1.3), or it could come from an unresolved collection of astrophysical 
sources. The cascade channel is useful for diffuse searches because, at the energies of 
interest for IceCube, the atmospheric electron neutrino flux is about a factor of 20 below 
the atmospheric muon neutrino flux (see section 1.2). In a search for a harder, diffuse source 
(for instance an E~ 2 extraterrestrial flux) on top of the atmospheric i?~ 3 - 7 neutrino flux, 
the new component should manifest itself as a break in the energy spectrum. Because the 
atmospheric background is lower for cascades, this break should show up at a lower energy 
in the cascade channel than the muon channel. With steeply falling fluxes, there should be 
more events at this lower energy and it should therefore be easier to resolve the spectral 
break. The better intrinsic energy resolution of cascades also makes the search for such 
a spectral break easier. Cascades deposit all of their energy in the detector and have an 
energy resolution in log(-E) ~ 0.2. Muons, on the other hand, deposit energy outside of the 
detector and have inherently worse energy resolution. 
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Chapter 3 

The IceCube Detector 

This chapter discusses the design and implementation of the IceCube detector, 
which uses the glacial icesheet at the South Pole as both a neutrino target and an optically 
clear medium in which to observe Cherenkov radiation. First, we discuss the optical proper- 
ties of glacial ice at the South Pole. An understanding of these optical properties is essential 
for event rate predictions from atmospheric and astrophysical neutrino fluxes and for the 
reconstruction of event characteristics like direction and energy. Next, we discuss the Ice- 
Cube array itself and the Digital Optical Module (DOM), IceCube's basic sensor element. 
Finally, we conclude with a description of drilling and sensor installation (deployment) and 
an update on the current status of construction. 

3.1 Optical Properties of South Pole Glacial Ice 

3.1.1 The Icesheet 

Unlike most particle physics experiments, IceCube does not make use of a well- 
controlled, man-made detector medium. In order to achieve the scale needed to observe 
extremely low astrophysical neutrino fluxes, IceCube has had to make use of a natural 
material — the glacial icesheet at the South Pole. This medium, having been subjected to 
the whims of earth's climatological past, is a complicated and rich one in which to embed 
a detector. 

At the South Pole, the icesheet is 2,820 m thick and was created by snow accu- 
mulation over a period of some 165,000 years [71]. Icesheets like this have been referred 
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to as "two-mile time machines" [72]. As one goes down in depth from the surface, one 
encounters layers of ice that were deposited as snow during long-distant periods of earth's 
climatological past. At the South Pole, a depth of 800-1,000 m corresponds to the onset of 
the Holocene, earth's current geological epoch which began ~1 1,000 years ago. At ~ 1,300 m 
depth lies the Last Glacial Maximum, the time period when the glaciers and icesheets last 
reached their maximum extent on earth ~20,000 years ago. The lowest deployed optical 
sensors of the array sit in ice which is ~90,000 years old [71]. 

At the very top of the icesheet lies the urn, a ~100 m layer of loosely packed snow 
and air. The firn extends until "pore closeoff", where the air no longer exchanges with the 
atmosphere. As overburden pressure increases with depth, the firn gives way to compactified 
poly crystalline ice. This ice contains impurities which were trapped in the snow when it 
fell, and the amount and type of impurity are related to the climatological conditions at 
the time of snowfall. For the purposes of IceCube, the most important impurity consists 
of micron-sized dust grains which were carried onto the icesheet as aeolian (wind-borne) 
aerosols. Many of these dust particles originated in earth's deserts and were transported 
around the upper atmosphere. The amount of this dust is related to how cold and windy a 
given climatological period is [73, 74] and varies strongly as a function of time, and therefore 
depth, in the icesheet. 

3.1.2 Basic Definitions 

Two optical properties of the icesheet are of particular importance for IceCube. 
The first is the absorption length A a , which is the distance over which a photon's probability 
of survival drops to 1/e. The second is the geometric scattering length X s , which is the 
average distance a photon travels between scatters. At each scatter, the photon has an 
angle change described by a scattering function p(8) which depends on the physics of the 
scattering interaction. The average of this scattering function is defined to be {cos 9). If 
(cos 9) is close to one, the scattering is highly peaked towards the forward direction. 

To incorporate the fact that many forward scatters are essentially equivalent to 
no scattering at all, we define an effective scattering length A e , which is the distance over 
which a photon has its direction randomized. To calculate A e , we make use of the following 
argument from reference [75] . A photon travels a distance \ s and is then scattered through 
an average angle (cos 9). It then travels another X s and is scattered again through an 
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average angle (cos 6). Up to this second scatter, the photon distance projected along the 
original direction of travel is A s + A s (cos#). Repeating this procedure over n scatters, the 
effective length of light travel projected along the original direction of travel is 

n 
i=0 

In the limit of large n, this becomes 

A - A * 

e 1-(COS0> 

If the scattering is forward-backward symmetric, (cos#) = and A e = X s . That is, in one 
geometric scattering length the direction is effectively randomized. For forward-scattering, 
(cos 6) > and so A e > A s . The result of the forward-scattering is that one has to travel an 
effective distance longer than the geometrical scattering length in order to fully randomize 
the direction. Calculations of Mie scattering (since the particle size is comparable to the 
wavelength of light) off micron-sized dust yield a value of {cos 9} ~ 0.94 [75]. 

It is often convenient to work in terms of the absorptivity a and the scattering 
coefficient b e rather than the absorption and scattering lengths themselves. These are 
defined by 



1 

a; 
i 

: aT 



Using pulsed and continuous light sources, the AMANDA collaboration measured the ab- 
sorption and scattering as a function of depth at the South Pole. Measurements of photon 
arrival time distributions were fit to full Monte Carlo simulations of photon propagation. 
The definitive discussion of the resulting optical properties map can be found in [75]. 

3.1.3 Scattering 

A full discussion of scattering mechanisms in glacial ice can be found in [76]. The 
two most important sources of scattering in South Pole glacial ice are submillimeter-sized 
air bubbles and micron-sized dust grains. 
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At depths shallower than ~1,400 m, submillimeter-sized air bubbles are trapped in 
the ice and are highly scattering. The mean free path monotonically increases with depth 
as the hydrostatic pressure compresses the bubbles. Eventually, as depth and pressure 
continue to increase, the air bubbles become unstable against a phase transition to a solid 
air hydrate phase. In this phase, which physically resembles ice, the gas molecules are 
trapped within the crystalline ice. By ~ 1,400 m depth at the South Pole, the fraction of 
bubbles transformed into clathrates has grown to 100% [77]. The refractive index of the 
clathrates almost perfectly matches that of the ice, making them essentially invisible. 

At depths where the air bubbles have disappeared, scattering becomes dominated 
by micron-sized dust impurities trapped in the ice. The main components of this dust 
are mineral grains, sea-salt crystals, liquid acid drops, and soot. Calculations suggest that 
scattering off dust should have a power law dependence on wavelength: 

b e oc \~ a 

where alpha is close to 1. Figure 3.1 shows the measured wavelength dependence of the 
scattering coefficient. 

Figure 3.2 shows the measured scattering coefficient as a function of depth below 
the surface of the icesheet at the South Pole for several wavelengths. At the depths of 
interest for IceCube, the effective scattering length varies from around 10 to 50 m. The 
dusty peaks labeled A, B, C, and D are related to known climatological features. 

3.1.4 Absorption 

A full discussion of absorption mechanisms in glacial ice can be found in [78]. For 
near-visible wavelengths, three absorption mechanisms are important for glacial ice at the 
South Pole. The absorptivity can be parameterized as [75] 

a(A) = A v e~ B v x + C7 dust A- K + A lR e' x ^ x 

The first and last terms are from intrinsic absorption processes in pure ice itself, and the 
middle term is due to absorption on dust grains. The first term is the so-called Urbach tail, 
which governs the short wavelength behavior of the absorption at wavelengths just above 
the electronic band-gap energy of the ice crystal. The last term gives an exponential rise in 
the red and infrared and comes from molecular stretching and bending modes of the pure 
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Figure 3.1: Wavelength dependence of the scattering coefficient b e . From [75]. 
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Figure 3.2: Scattering coefficient 6 e as a function of depth below the surface of the icesheet. 
From [75]. 
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Figure 3.3: Absorption coefficient a as a function of depth below the surface of the icesheet. 
From [75]. 



ice crystal. The power law exponent of the second term due to dust was measured to be 
k = 1.08 ±0.01. 

Figure 3.3 shows the measured absorptivity as a function of depth below the surface 
of the icesheet at the South Pole for several wavelengths. At the depths of interest for 
IceCube, the absorption length is around 75 to 100 m. At a wavelength of 532 nm, the 
intrinsic absorption of ice dominates the absorption from dust. The linear rise with depth 
at this wavelength reflects an increase of ice temperature with depth and the corresponding 
increase of the intrinsic absorption with temperature. At the other wavelengths the same 
climate features as in the scattering coefficient can be observed. 

Figure 3.4 shows the measured wavelength dependence of the absorptivity at sev- 
eral depths. Note that the glacial ice at the South Pole is an order of magnitude cleaner 
and more transparent between 200 nm and 400 nm than ice made in a laboratory. 

Figure 3.5 shows the full three-dimensional dependence of scattering and absorp- 
tion function of depth and wavelength. 
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Figure 3.4: Wavelength dependence of the absorptivity a for glacial and laboratory ice. 
From [75]. 




Figure 3.5: Scattering and absorption coefficients as a function of depth and wavelength. 
From [75]. 
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3.1.5 Diffusive Regime 

In the early days of AMANDA, when sensors were deployed in very bubbly ice 
with scattering lengths on the order of 10 cm, it was realized that photon transport could 
be described by a diffusion equation with absorption. For an isotropic point source of light, 
the photon arrival time distribution at a distance d away is given by [76] 

n(d ' t] = jidw exp (5^) exp (^) 

where the diffusion constant D is given by D = ^ If we integrate this equation over 
time to get the total number of photons expected at a distance d from the source, we have 

N(d) oc \(T dlXv 

where the propagation length X p = \J Xe ^ a is typically around 25 m. This approximation 
is generally valid in the regime where d > 5A e and has been used in a cascade energy 
reconstruction which will be discussed in the next chapter. 

3.2 Hardware 

In essence, the IceCube detector is simply an array of photosensors embedded in 
the deep glacial ice at the South Pole. These photosensors detect the Cherenkov light emit- 
ted by charged secondary particles produced in high energy neutrino interactions with water 
molecules. In practice, the hardware needed to accomplish this task, both in the ice and on 
the surface, is quite complicated and has required significant research and development. The 
digital technology and robust, low-power, low-maintenance design of the IceCube experi- 
ment is the culmination of the experience gained by the AMANDA collaboration, which 
spent more than a decade building and operating a prototype neutrino telescope at the 
South Pole. 

IceCube's basic sensor element is the Digital Optical Module (DOM), a photo- 
multiplier tube and its accompanying electronics housed in a glass pressure sphere. DOM's 
are connected to their neighbors and to the IceCube Lab at the surface via a bundle of 
twisted-pair copper wires. A "string" consists of 60 DOM's and their cable bundle and is 
deployed in a vertically-drilled hole in the glacial ice. DOM's are separated from each other 
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by 17 m on each string, and the string itself is lowered so that all of the DOM's lie at depths 
between 1,450 m and 2,450 m where the optical properties of the ice are most favorable. 

IceCube's baseline plan calls for 80 strings with a nominal inter-string spacing of 
125 m. The 17 m DOM spacing and 125 m string spacing leads to a neutrino energy thresh- 
old ~100 GeV. Recently, approval was granted to install 6 additional densely-instrumented 
strings at the very center of the array to form an inner detector called DeepCore [79] . Each 
DeepCore string will host 50 DOM's separated from each other by 7 m at depths below 
2,100 m where the ice is the clearest. 

In addition to this "in-ice" detector, IceCube also has a surface array for studying 
cosmic ray air showers. IceTop is made up of pairs of frozen water tanks located around 
25 m from the top of each in-ice string. The two tanks are separated from each other 
by 10 m, and each tank has two DOM's for sensing Cherenkov light from the passage of 
relativistic charged particles from air showers. 

Figure 3.6 shows a schematic of the IceCube detector. 
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Figure 3.6: Schematic of the IceCube detector. 
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Figure 3.7: The Digital Optical Module (DOM). 

3.2.1 The Digital Optical Module (DOM) 

Reference [80] gives a detailed description of the hardware and electronics found 
in the IceCube Digital Optical Module (DOM). A DOM schematic and photograph are 
shown in figure 3.7. The DOM's light-sensing element is a 10 in. photomultiplier tube 
(Hamamatsu R7081-02) that achieves a gain of 10 7 over its 10 dynode stages running at a 
supply voltage around 1500 V. Its bialkali photocathode is sensitive between 300 nm and 
600 nm with maximum sensitivity around 420 nm. Figure 3.8 shows the measured PMT 
wavelength acceptance (the product of the quantum efficiency and the collection efficiency) 
as a function of wavelength. Dark noise rates are typically around 540 Hz. 

The PMT and its accompanying electronics are housed in a 13 in. diameter, 0.5 in. 
thick glass pressure sphere so that they can withstand the high pressures (up to 10,000 PSI) 
in the water column during deployment and refreezing. The glass was formulated specifically 
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Figure 3.8: Wavelength acceptance of the PMT (top) and glass sphere (bottom). 
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for IceCube to maximize the transmission in the region of PMT sensitivity. The measured 
glass transmission as a function of wavelength is also shown in figure 3.8. The PMT is 
mechanically and optically coupled to the glass with a silicone gel (GE6156). 

The current pulse from the PMT travels down several parallel paths. One part is 
sent to a discriminator which decides whether or not the DOM has triggered. The other 
part is sent down a 75 ns delay line while the discriminator makes it choice. Typically the 
discriminator threshold is set to 0.25 photoelectrons. If the discriminator fires, it initiates 
digitization of the signal by several analog-to-digital converters at the end of the delay line. 

The Analog Transient Waveform Digitizer (ATWD) [81] is a low-power, fast- 
sampling ADC. It samples at 300 MSPS over a time window of 400 ns. Three different 
amplifiers are placed before the ATWD inputs (xl6, x2, and x0.25) to provide large dy- 
namic range. If any sample in the xl6 channel saturates, the x2 channel is digitized. If any 
sample in that channel saturates, then the x0.25 channel is digitized. These channels are 
referred to as ATWD0, ATWD1, and ATWD2. Offline analysis combines these channels 
into one calibrated waveform. Each DOM has two separate ATWD's, and it ping-pongs 
back and forth between them to minimize dead time. 

The second ADC, the flash ADC (FADC) provides sparser sampling over a longer 
time window. It is meant for particularly long, bright signals. It samples at 40 MSPS over 
6.4 fxs. 

Each DOM is equipped with 12 on-board "flasher" LED's that operate at 405 nm. 
The flashers are spaced around the azimuth of the DOM, 6 pointing in the horizontal 
direction and 6 tilted upwards at an angle ~ 45°. Both the total brightness and the length 
of the flash are configurable. The flashers are used for calibration and verification tasks 
(measuring ice properties, verifying the timing and geometry of the array) as well as more 
physics-centric tasks. Since the flashers are essentially point sources of light, they are useful 
for verifying the position and energy reconstruction algorithms that are used for neutrino- 
induced cascades. 

3.2.2 Local Coincidence 

Each DOM also has circuitry for implementing a local coincidence requirement 
between neighboring DOM's. When a DOM's discriminator fires, it sends signals directly 
to its two nearest neighbors to see if they have also fired. These DOM's can relay the signal 
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on to their neighbors to probe the behavior of more distant DOM's on the string. Typically, 
local coincidence conditions are imposed that require at least one of a DOM's two nearest 
neighbors or two next-to-nearest neighbors to have fired in order for data from that DOM 
to be sent to the surface. The local coincidence time window is typically ±l//s around the 
launch of a DOM. 

Local coincidence significantly reduces the number of accidental noise hits that 
enter the datastream. If a DOM fires at r = 540 Hz from dark noise, we can calculate 
the probability that one of its four nearest neighbors will also have a noise hit in the local 
coincidence window t = 2fis. The probability that there are no hits in any of the four 
DOM's is 

( e - ri ) 4 = e" 4ri « 1 - Art 
The probability of at least one hit in those four DOM's is thus 

1 - (1 - Art) = Art 

So the rate at which a DOM and at least one of its four nearest neighbors trigger from noise 
is Ar 2 t. In the case of the numbers given above this is ~ 0.6 Hz, three orders of magnitude 
below the dark noise rate. 

3.2.3 Triggering 

IceCube's triggering system is very flexible and can accommodate multiple triggers 
at the same time. It can be configured to use a simple multiplicity trigger (SMT), which fires 
when a given number of DOM's are hit within a given time window. It can be configured to 
use a so-called string trigger, which fires when a given number of DOM's on a single string 
are hit within a given time window. The string trigger is useful for low-energy atmospheric 
or WIMP-annihilation neutrinos. The triggering system also has the ability to search for 
more complicated event topologies in real-time. For the neutrino-induced cascade analysis 
in this dissertation, the simple majority trigger was used that required 8 hit DOM's in a 
time window of 5 /is. 
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3.3 Drilling and Deployment 

Each string is deployed in a vertically-drilled hole in the glacial ice. Drilling 
proceeds in two stages. First, a firn drill melts through the loosely-packed snow and air 
that make up the ~100 m thick firn. It operates by running heated fluid through a closed 
loop of coils which are in contact with the snow. Hot-water drilling, which is used for the 
deep ice drilling, is not appropriate for the firn as water leaks into the loosely packed snow 
rather than remaining in the hole. 

After firn drilling, the remainder of the hole is drilled by IceCube's dedicated 5 MW 
hot-water drill. An array of heaters is used to heat water, which is then pumped to the 
drill-head. The drill-head is attached to a heavy weight stack and shoots a high-pressure 
jet of hot water into the bottom of the hole. A pump returns water back to the heaters for 
re-heating. 

The entire drilling process takes 2 to 3 days. When it is complete, a column of 
water is left in the hole, which is approximately half a meter in diameter. Next, a team of 
deployers makes the mechanical and electrical connections for each of the 60 DOM's to the 
string, and the string is then lowered into the hole. Deployment takes ~ 12 hours and is 
completed well before the hole begins to refreeze. The final depth of the string is measured 
by several pressure sensors which measure pressure in the water column. 

Figure 3.9 shows the 5 MW hot-water drill in action. In the lower left of the image, 
a series of red trailers houses the heaters and pumps. Hot water flows from these trailers 
through the black hose which snakes up to the tower structure in the upper center of the 
image. This tower is situated over the hole which was being drilled at the time. 

3.4 Status 

During the austral summer of 2004-2005, the first IceCube string was successfully 
deployed and commissioned at the South Pole. Since then, the pace of deployment has been 
rapidly accelerating with each season. In the 2005-2006 season, 8 strings were installed and 
in 2006-2007, 13 more were deployed for a total of 22. This 22-string array, known as IC- 
22, is the array used for the analysis in this dissertation. The 2007-2008 season saw the 
installation of 18 more strings, bringing the array up to 40 strings. IC-40 just completed 
its one year data-taking phase. This past season, 2008-2009, 19 strings were deployed. 



54 




Figure 3.9: IceCube's 5 MW hot-water drill at work. In the lower left, a series of red trailers 
houses heaters and pumps. Hot water flows from these trailers through the black hose which 
snakes up to the tower structure in the upper center. This tower is situated over the hole 
which was being drilled at the time. 
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Currently, IceCube is taking data in its IC-59 configuration. The full array is due for 
completion in 2010-2011. 

Figure 3.10 shows a surface map of the string layout of IC-22. Several salient 
features of this string layout are important for a search for neutrino-induced cascades. As 
will become clear later, the main background for a neutrino-induced cascade search consists 
of muons which have a large radiative energy loss inside the detector that mimics a cascade. 
In order to veto these events, it is desirable to have a contained region of the detector 
which is surrounded by many strings to capture early light from the muon. For IC-22, 
there are only 6 strings (57, 48, 39, 66, 58, and 49) which have a single layer of strings 
surrounding them. This means that the veto region is small, and vetoing muon events with 
large radiative losses will be a challenge. Second, the geometry is quite irregular, with 
several strings (78, 72, and 46) jutting out from the main body of array. These strings are 
likely to be less useful for a cascade search, and so the effective size of the detector will be 
reduced. 
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Figure 3.10: Surface map of the 22-string IceCube array (IC-22). 
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Chapter 4 



Reconstruction of 



Neutrino-Induced 



Cascades 



As we have seen, a particle shower or 



cascade" from a charged-current or neutral- 



current neutrino interaction is a relatively local deposition of energy. In a sparsely instru- 
mented detector like IceCube, these cascades look effectively like point sources of Cherenkov 
light. For each candidate cascade event in our detector, we would like to reconstruct its 
physical parameters. This chapter discusses the methods for doing so. Event reconstruction 
will allow us to separate neutrino-induced events from air shower muon background events 
and to make sense of the resulting neutrino sample. 

4.1 Reconstruction Basics 

A neutrino- induced cascade is characterized by a set of unknown parameters {a}: 



where f v is the cascade vertex position, t v is the vertex time, E is the visible energy of the 
cascade, and h is its direction. To reconstruct these parameters from our observational data, 
we make use of maximum likelihood estimation [29] . For a given set of observational data 
{x}, we ask which cascade parameters {a} maximize a given likelihood function Jzf(x|a). 
If the Xi components of the data {x} are independent, we can write 



{a} = {r v ,t v ,E,h} 
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Jf(x|a) = JJp(xj|a) 

i 

where p(xj|a) gives the probability of our experiment measuring the quantity Xi given the 
cascade parameters {a}. For example, p(xi\a) could tell us the probability of observing 
a single photoelectron at a given time in a single DOM from a cascade with parameters 
{a}. Or, it could tell us the probability of observing a total of N photoelectrons in a single 
DOM from a cascade with parameters {a}. Maximizing the likelihood function Jzf(x|a) 
corresponds to finding the set of cascade parameters that makes the observed data most 
likely to have occurred. 

Various likelihood functions can be constructed for neutrino-induced cascades. 
The most important uses photoelectron timing information to reconstruct the position and 
time of the cascade vertex. This is discussed in section 4.3.1. In addition, two different 
likelihood functions are available to reconstruct the visible energy of the cascade. The 
traditional method does not take into account the depth-dependent optical properties of 
the ice because of the conceptual and technical difficulties of doing so. It is described in 
section 4.3.2. For this analysis, a better performing depth-dependent energy reconstruction 
was developed which is described in section 4.4. 

Directional reconstruction of cascades has proven a difficult challenge. While the 
light emission from a cascade starts out with a directional anisotropy due to the Cherenkov 
light emission, scattering in the ice tends to isotropize it very quickly. For this reason, 
directional reconstruction was not performed at any stage of this analysis. Very recent 
work on improved reconstruction methods suggests that it may be possible to achieve a 
directional resolution of 30°-35° for 10 TeV-10 PeV cascades in the near future [82]. 

In practice, it is numerically easier to minimize a function rather than to maximize 
it. Consequently, we work with the negative log of the likelihood function instead of the 
the likelihood function itself. In addition, because the likelihood landscapes are often quite 
complex, it is desirable to have a set of seed parameters {a} that we believe lies in the 
vicinity of the true global minimum of the negative log likelihood function. This seed is 
usually provided by a quick-running "first-guess" algorithm. The first-guess algorithm for 
the cascade vertex and time which is used to seed the likelihood reconstructions is described 
in the next section. 
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4.2 The cfirst First-Guess Algorithm 

The cfirst program provides a first guess for the cascade vertex by making an 
analogy between photoelectrons observed in DOM's (also known as "hits") and masses 
distributed in space [29]. We take as our first guess the center-of-gravity (COG) of the hits 
defined by: 

cog ee ^>*r • n 

1=1 

where a is the number of photoelectrons recorded by a DOM at position fl , w is an arbitrary 
amplitude weighting factor, and iV c h is the number of "hit" DOM's (that is, the number of 
DOM's which received photoelectrons in the event). 

In general, this first guess works well. In the x and y directions the COG is 
symmetric about the true vertex position with a resolution of around 25 m. In the z 
direction, the 17 m resolution reflects the denser DOM spacing along the string. However, 
in the z direction the COG tends to lie 10 m above the true vertex. This is because a 
DOM's PMT points downwards, and so the DOM's above a cascade see more light than the 
DOM's below . The z value of the first guess also tends to be poor near layers of dirty ice 
in the detector, since in these regions light may be preferentially absorbed in one direction. 
Finally, the COG must by definition lie within the instrumented volume, so it is not a good 
guess for events which are truly outside of the detector. 

After the COG position is found, it is used to calculate a first guess for the vertex 
time. The time calculation relies on the fact that hits very close to the cascade vertex have 
a good chance of having arrived unscattered or close to unscattered. We choose our vertex 
time to maximize the number of these direct hits. 

For a hit at time ti at a distance di from the COG, we start by assuming that it 
came directly from the COG and therefore take as our trial vertex time i v = ti — Next, 
we define a time residual variable which tells us how delayed from direct travel another hit 
tj is, given our trial vertex time: 

tr = tj ^direct = tj (t v H — ) = tj (ti I ) = (tj ) (t{ ) 

Cice Cice Qce C; ce C; ce 

Positive values of t r indicate that the light was scattered on its way from the COG to the 
hit DOM and so arrived later. We calculate the residuals for all of the other hits in the 
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event given our trial vertex time and count how many of them lie in a [0 ns, 200 ns] window. 
These are the hits that are close to direct for our trial vertex time. 

We repeat this procedure taking the trial vertex time from each hit in the event. 
Our final vertex time is chosen to be the earliest trial time that results in at least iVtrigg = 3 
hits in our [0 ns, 200 ns] direct time window from this procedure. In case the N tr ig g = 3 
condition is not met, the vertex time is estimated as the time of the first hit in the event. 
This first guess method for the vertex time has a resolution of around 90 ns, though by 
ignoring scattering it tends to push the vertex time as 35 ns too late. 

4.3 Log-Likelihood Reconstructions 

This first-guess for the cascade vertex and time can then be used to seed a full 
likelihood reconstruction of the position, time, and energy. 

4.3.1 Vertex Reconstructions 

The vertex likelihood functions use photoelectron arrival time information to re- 
construct the cascade position and time. The individual pdf's p(xj|a) are parameterized in 
terms of the time residual defined earlier: 

tr = tj ^direct 

The pdf p(t r \a) gives the probability of measuring a single photoelectron with time residual 
t r in a given DOM from a cascade with parameters {a}. 

If South Pole glacial ice were not a strongly scattering medium, this time residual 
pdf would be a simple delta function at t r = 0. However, scattering and absorption in the 
ice significantly complicate the situation. For a vertex at a given position in the detector, 
light will propagate through a layered system of scatterers and absorbers. Very close to 
the vertex, we still expect the residual function to peak around zero, as there is still a high 
probability of unscattered light. But now the distribution should have a long tail towards 
large residuals as well. Far from the vertex, we expect very little unscattered light. 

Traditionally, AMANDA and IceCube researchers have used an analytic function 
called the "Pandel function" as a convenient parameterization of the time residual pdf as a 
function of the distance d from the cascade vertex [83]. The Pandel function is given by: 
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Figure 4.1: Pandel time residual pdf (dashed) compared to photonics [84] Monte Carlo 
simulation (solid) at distances of 25 m (left) and 125 m (right) from a cascade vertex located 
at a depth of 1948 m. 
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where the normalization N(d) is 

-d/X 



N(d) = e~ d/Xa ^1 + 



rc-u 



Here, A a is the absorption length, A and r are free parameters, and T(x) is the standard 
gamma function. The free parameters A and r are determined by fitting this functional 
form to time residual distributions from a full Monte Carlo simulation of scattering and 
absorption in layered glacial ice. The values used for the reconstructions in this analysis 
are r = 450 ns, A = 47.0 m, and A a = 98.0 m. They have been chosen to fit the widest 
range of time residual distributions and are not related to actual physical parameters. 

Time residual distributions at two different distances are shown in figure 4.1. For 
small distances, the factor t r w^~ causes the probability to rise sharply to a pole at zero 
time residual. This reflects the high probability for unscattered photons at short distances 
from the vertex. At a critical distance d cr ;t = A, the exponent of this factor changes sign 
and the distribution no longer peaks at zero time residual. This reflects the fact that there 
is very little chance of getting direct, unscattered light far from the vertex. 

Two additional features must be added to this pdf to make it a full description of 
the detector. First, the DOM PMT and electronics don't have perfect timing resolution. 
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Table 4.1: Vertex position resolutions in meters at various levels of this analysis weighted 
to an atmospheric energy spectrum. 



Analysis Level 


X 


y 


z 


Level 3a 


20.87 


19.17 


10.30 


Level 4a 


11.48 


10.84 


6.75 


Final Level 1 


10.02 


9.85 


5.09 


Final Level 2 


9.71 


9.58 


5.47 


Final Level 3 


9.37 


9.29 


5.82 



To account for this, we smear the Pandel residual distribution above with a gaussian rep- 
resenting our finite timing resolution. The gaussian was chosen to have mean /x = 0ns and 
a = 15 ns. Second, we account for noise hits by adding a small, flat noise probability term 
of 10~ 10 to the pdf. 

The final pdf incorporating all of these elements is known as the UPandel function 
Pu(t r \d) and has the following form: 



N(0,a) fort r <0 

pu(t r \d) = { P 3 (t r \d) for < U < y/2Tro 

p(t r \d) for t r > \p2^o 

where iV(0, a) is the gaussian representing our timing resolution, p(t r \d) is the normal Pandel 
function defined above, and Ps(t r \d) is a third-order polynomial interpolation function. Our 
full likelihood function is then: 



^UPandel (x| a) = JJ £>[/ (t r | a) 
hits 

Two versions of the UPandel vertex reconstruction were run as part of this analysis. The first 
used only the earliest hit recorded in each DOM. A more time-intensive reconstruction used 
all of the hits recorded in each DOM and was the best-performing fit. Position resolutions 
from the all-hit reconstruction are shown in figures 4.2-4.4 for various levels of this analysis. 
See chapter 6 for a definition of the analysis levels. Table 4.1 lists the vertex resolutions as 
a function of analysis level. 
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Figure 4.2: x vertex position resolutions for E 2 (dashed) and atmospheric (solid) spectra 
at level 3a (top), level 4a (middle) and final level 1 (bottom) of this analysis. 
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Figure 4.3: y vertex position resolutions for E 2 (dashed) and atmospheric (solid) spectra 
at level 3a (top), level 4a (middle) and final level 1 (bottom) of this analysis. 
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Figure 4.4: z vertex position resolutions for E 2 (dashed) and atmospheric (solid) spectra 
at level 3a (top), level 4a (middle) and final level 1 (bottom) of this analysis. 
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4.3.2 PHit-PNoHit Energy Reconstruction 

Because it relies on photoelectron arrival times, our UPandel vertex reconstruction 
only takes advantage of information from DOM's which were hit. However, there is also 
information contained in the fact that a DOM was not hit. This information is incorporated 
into the "PHit-PNoHit" energy reconstruction. 

To construct the PHit-PNoHit likelihood function, we need to know the probability 
that a DOM is hit given a cascade vertex and energy. To do so, we make several assumptions. 
First, we assume that the cascade is an isotropic point source of light. Next, we make the 
diffusive scattering approximation (see section 3.1.5). Finally, we make use of the facts 
that the total track length and therefore the total light intensity are proportional to the 
cascade energy. Given these assumptions we can write an expression for the mean number 
of photoelectrons at a distance d from the cascade [58]: 

a 

where Iq is a normalization constant determined from fits to Monte Carlo simulations and 
the propagation length is defined by X p = yj Typically, Iq = 1.4 GeV^ 1 m and 
X p = 25 m. 

From Poisson statistics, we know that the probability of observing zero photoelec- 
trons is e~^. Therefore we can write: 



Pease -i pease -i &~t J > 

hit — 1 nohit — 1 " 

This hit probability for several cascade energies is plotted in figure 4.5. 

To incorporate noise, we add a constant noise hit probability term: 



p -i p pease , p 

-'hit — 1 -'nohit — -» hit t ■'noise r hit -'"noise 

where the subtracted term avoids double counting the source of a hit. Using this expression 
for Phit) we have our final PHit-PNoHit energy likelihood function: 



■^PHit-PNoHit 

(E)= Put(E,d)x fUitOM) 

all hit DOM's all un-hit DOM's 



67 




distance [m] 

Figure 4.5: From left to right, hit probability for a 100 GeV, 1 TeV, 10 TeV, and 100 TeV 
cascade as a function of distance from the cascade vertex. 

It is traditional in neutrino astronomy to express energy resolution in terms of 
log(£J reco ) — log (Eligible)- In AMANDA, the PHit-PNoHit reconstruction resulted in an 
energy resolution between 0.1 and 0.2 in log(E). 

This energy reconstruction method has several drawbacks, however. First, it as- 
sumes that we are in the diffusive regime, which may not hold for cascades relatively close 
to a DOM or a string. Second, by assuming an average propagation length, it ignores the 
dramatic depth variation of the optical properties of the ice. A cascade at one depth can 
result in significantly more or less detected light depending on the local ice properties. For 
example, a parallel IC-22 cascade analysis which made use of the PHit-PNoHit energy re- 
construction found that cascades at the top of the detector were shifted in energy by -0.5 
in log(E) while cascades at the bottom of the detector were shifted in energy by +0.5 in 
log(E). 

To overcome these issues, we developed a depth-dependent energy reconstruction 
discussed in the next section. 
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4.4 Depth-Dependent Energy Reconstruction 

In conjunction with a colleague, I developed and implemented a fast, analytic cas- 
cade energy reconstruction that takes the depth variation of the optical properties of glacial 
ice into account. This reconstruction program, AtmCscdEnergyReco, takes a reconstructed 
cascade vertex and outputs a reconstructed energy given that vertex hypothesis. 

To account for the depth variation of the optical properties, AtmCscdEnergyReco 
makes use of photon tables generated by a dedicated Monte Carlo package called photonics 
[84]. photonics tracks photons from the cascade through the dust layers of the glacial ice, 
simulating the full wavelength dependent scattering and absorption along the way. The 
result is a set of tables which give photon amplitudes and time delay distributions at a 
given location in the ice from a given cascade vertex position. These tables have built into 
them assumptions about the refrozen ice which surrounds the IceCube DOM's as well as 
the glass transmission and PMT quantum efficiency. To get the full amplitude, the table 
amplitude is scaled linearly by the cascade energy. 

These tables, known as "photorec" tables, tell us the expected number of photo- 
electrons in a given DOM for a cascade located at a particular depth in the detector. They 
are used in the likelihood functions described below. 

4.4.1 AtmCscdEnergyReco: Noiseless 

If we ignore noise hits, we can build an energy likelihood function as a product of 
Poissonian terms for the total charge seen by a DOM: 

* {S) = n ME) " P T" (E ' 

Hit+UnHit DOM's pc ' 

where fJ>(E) is the expected number of photoelectrons in a given DOM from a cascade of 
energy E and n pc is the observed number of photoelectrons. Taking the negative log of this 
expression we have: 

-logJSf (E) = - J2 (npelog(M(£0) " K E ) ~ log(npeO) 
Differentiating with respect to E and setting the result equal to we have: 
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Remembering that the amplitude scales linearly with the cascade energy E, we can write 
fji(E) = no(f v , tdom)E where /J>o(r v ,rDOM) is the amplitude from the photorec table which 
depends on the cascade vertex position r v and DOM position tdom but is independent 
of energy. E is the cascade energy in GeV. Using this expression we can write the above 
equation as 

Y] Ho(r v ,r DO M) = y~] — 7=; ^ x fJ>o(r v ,r DO M) 

^ ^ Ho{r v ,r D OM)E 

which implies that 



£[GeV] 



YlfJ>o(r v ,r DO M) 

So to reconstruct the cascade energy, we simply sum the observed charges and divide by 
the sum of the photorec table amplitudes for the given cascade vertex hypothesis. 
Finally we can calculate a reduced negative log-likelihood 

rlogJzr = 



iVch 

where _2? m i n is the likelihood function evaluated at the reconstructed energy. For a "perfect" 
cascade where each DOM sees exactly as much charge as the photorec tables say that it 
should, the reduced negative log-likelihood tends to one. This parameter should be useful 
for separating signal from background as muons with extra hits should have worse (that is, 
larger) reduced negative log-likelihood parameters. 

4.4.2 AtmCscdEnergyReco: With Noise 

With the addition of a noise term, the expected charge becomes [i{E) + rT where 
r is the noise rate in Hertz and T is the readout window length in seconds. Following the 
same logic as above, we arrive at the following equation: 
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~ 2 #T E ^t + Yl fboAf) = 

i !M>{r V! pDOM) i 

The presence of the noise term means that the solution for the energy is no longer ana- 
lytic. However, it can be found by an iterative root finding algorithm. A unique root and 
hence a unique solution for the energy always exists. The root-finding algorithm has been 
implemented with the GSL library [85]. 

4.4.3 Performance 

AtmCscdEnergyReco is currently the best-performing cascade energy reconstruc- 
tion in use by the IceCube collaboration. Energy resolution plots for various levels of this 
analysis are shown in figure 4.6. Figure 4.7 shows a scatter plot of reconstructed energy as 
a function of visible energy in the ice. Table 4.2 lists the energy resolutions as a function 
of analysis level. 

Because it makes use of the photorec tables which have knowledge of the depth 
dependence of the optical properties of the glacial ice, AtmCscdEnergyReco performs well 
as a function of depth in the detector. Figure 4.8 compares the energy resolution in the 
top, middle, and bottom of the detector. 

Finally, studies were carried out to verify the performance of AtmCscdEnergyReco 
using the in-situ flasher LED's located in each DOM. Figure 4.9 shows the reconstructed 
number of photons per flash for LED's located at various depths in the ice for this recon- 
struction and the PHit-PNoHit reconstruction. Our AtmCscdEnergyReco algorithm does a 
much better job at reducing the influence of local optical properties on the reconstructed 
energy. 
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Table 4.2: Energy resolutions in log(E) at various levels of this analysis weighted to an 
atmospheric energy spectrum. 



Analysis Level 


RMS of log(£ rcco ) - log(£ visiblc ) 


Level 3a 


0.26 


Level 4a 


0.23 


Final Level 1 


0.16 


Final Level 2 


0.17 


Final Level 3 


0.17 
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Figure 4.6: From top to bottom, energy resolutions for E 2 (dashed) and atmospheric 
(solid) spectra at level 3a, level 4a, and final level 1 of this analysis. 
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Reconstructed Vs. True Energy at Level 4a 




Figure 4.7: Reconstructed energy vs. visible energy in the ice at level 4a of this analysis. 
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Energy Resolution at Level 4a (Atmo. Spectrum) 
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Figure 4.8: Energy resolution at the top (red), middle (blue), and bottom (green) of the 
detector at level 4a of this analysis for an atmospheric spectrum (top) and an E~ 2 spectrum 
(bottom). 
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Figure 4.9: Reconstructed number of photons per flash for flasher LED's located at vari- 
ous depths in the ice from [86] for AtmCscdEnergyReco (top) and PHit-PNoHit (bottom). 
The broken line shows the true number of photons per flash. AtmCscdEnergyReco does a 
much better job of smoothing out the effects of local ice properties on the reconstructed 
energy. The residual structure may indicate remaining deficiencies in the modeling of the 
ice properties. 
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Chapter 5 

Simulation 

Analysis in IceCube is very heavily dependent on Monte Carlo simulations to 
model both the muon background from cosmic ray air showers and the neutrino signal. This 
chapter discusses the simulation of both classes of events. First, air shower and neutrino 
generation will be described. Then, the detector simulation chain, common to both types 
of simulation, will be explained. Figure 5.1 is a schematic of the various programs that go 
into simulating signal and background events. 

5.1 Background Simulation 
5.1.1 CORSIKA 

The main background for all IceCube analyses consists of downgoing muons from 
cosmic ray air showers. To simulate this background, we utilize the Monte Carlo simulation 
program CORSIKA (COsmic Ray Simulations for KAscade, version 6720) [87]. 

CORSIKA propagates cosmic ray primaries through the atmosphere to their point of 
interaction with an air nucleus. Various high energy hadronic interaction models can then 
be used to simulate the interaction, including VENUS [88], QGSJET [89], DPMJET [90], 
or SIBYLL [91]. For this analysis, SIBYLL was used. The secondary particles produced in 
the interaction are propagated and allowed to further interact or decay. 

The background simulations utilize a primary cosmic ray energy spectrum known 
as the Hoerandel polygonato model (Greek for "many knees") [92]. In the polygonato 
model, the energy spectrum of each primary element (p, He, C, N, O, etc.) consists of two 
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Figure 5.1: Schematic of the various programs that go into simulating signal and background 
events. Event generators are shaded blue and the elements of the detector simulation are 
shaded green. 
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power laws patched together at a cutoff energy called the knee. The cutoff energy is an 
increasing function of nuclear charge Z. Individual primaries up to iron were simulated from 
600 GeV to 10 11 GeV and then added together according to the Hoerandel model. 

The simulation of each air shower results in a list of muons at the surface of the 
earth. These muons are then handed off to the detector simulation which is described in 
section 5.3. 

5.1.2 Energy Weighting 

The cosmic ray primary spectrum is a power law E"" 1 with index 7 ~ 2.7 steep- 
ening to 7 ~ 3 around an energy of 4 PeV. Experience from previous analyses as well as 
preliminary studies for this analysis have suggested that primaries from the higher energy 
end of the spectrum are more likely to result in the energetic muons with large radiative 
losses that survive to high cut levels of a cascade analysis. Consequently, it is more efficient 
to simulate a harder energy spectrum rather than the natural cosmic ray spectrum in order 
to oversample these higher energy air showers. 

The bulk of the Monte Carlo used for the analysis in this dissertation was generated 
with a primary energy spectrum that was one power harder than the natural spectrum. 
Event weights were stored and later applied to reweight back to the natural spectrum. 

5.1.3 Coincident Events 

In a detector as large as IceCube, there is a non-negligible probability that muons 
from two or more independent air showers may enter the detector volume together. These 
so-called coincident events are especially pernicious for muon neutrino analyses. Muons 
from one shower may cause hits at the bottom of the detector while muons from the other 
shower cause slightly later hits in the upper part, thereby mimicking an upward-going 
neutrino-induced muon. 

To find ways to eliminate this class of background, doubly coincident and triply 
coincident air shower events were simulated. First, individual CDRSIKA air showers were 
simulated, and events which resulted in at least one hit DOM in the detector were written 
to a file. Then, these events were paired together with a relative time difference less than 
5t = 21/is. A calculation was then carried out to assign the proper livetime to this collection 
of coincident events. That is, if we have N showers and we construct N/2 pairs, we need to 
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calculate the true amount of time we would have to wait to see N/2 coincident events [93]. 
5.1.4 Livetime 

CDRSIKA air shower simulation is an extremely computationally intensive task. 
Even with the resources of the entire IceCube collaboration running for ~ 6 months, we were 
only able to generate ~ 20 days of effective detector livetime. Various optimization schemes 
were tried but ultimately proved unsuccessful on the timescale of this atmospheric neutrino- 
induced cascade analysis. Producing enough background Monte Carlo will continue to be a 
significant challenge for future cascade analyses and is currently a very active research area 
within the IceCube collaboration. 

Appendix A has details on the calculation of effective livetime for a set of weighted 
CORSIKA events. 

5.2 Signal Neutrino Simulation 

5.2.1 neutrino-generator 

Neutrino simulation is done within the IceCube simulation framework IceSim and 
is handled by the neutrino-generator program, which in turn is based on the ANIS (All 
Neutrino Interaction Simulation) code [94]. neutrino-generator was used to generate 
neutrino events distributed uniformly over both hemispheres of the earth. The neutrinos 
were then propagated through the earth taking into account earth absorption (where the 
neutrino suffers a charged-current interaction somewhere deep inside the earth) and neutral- 
current regeneration (where the neutrino undergoes a neutral-current interaction resulting 
in a lower energy neutrino). For electron antineutrinos, the resonant scattering off electrons 
(the Glashow resonance) is also included. The Preliminary Reference Earth Model [95] 
is used to model the density profile of the earth and is shown in figure 5.2. The CTEQ5 
(Coordinated Theoretical-Experimental Project on QCD) parton distribution functions [54, 
55] were used for the neutrino cross sections which are also shown in figure 5.2. 

When the neutrino gets within a certain distance of the detector, an interaction 
of a given type is forced to occur, and the probability for that interaction is stored for 
event reweighting. In order to get absolute event rates for the atmospheric neutrino flux, 
the probability weight is combined with appropriate flux weights to take into account the 
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Figure 5.2: Earth density model and neutrino cross sections used by neutrino -generator. 
The cross section also includes the Glashow resonance, which is not shown in the figure. 

atmospheric flux variation as a function of energy and zenith angle. 

The output particles from the forced neutrino interaction are handed off to the 
CMC program, described in the next section. The output of CMC, in turn, is handed off to 
the detector simulation which is described in section 5.3. 

For this analysis, 4,930 electron neutrino files and 2,494 muon neutrino files were 
simulated. Each file contained 200,000 generated events from 40 GeV to 10 10 GeV with an 
E~ 2 energy spectrum for the primary neutrino. The interaction vertices were distributed 
around the detector in a cylinder of radius 1200 m and length 1700 m. 

5.2.2 CMC 

We often refer to cascades as point sources of light because of the sparse instru- 
mentation of IceCube. In reality, cascades do have a longitudinal profile and at energies 
above a TeV can extend over several meters. To simulate the longitudinal development of 
electromagnetic and hadronic cascades, the program CMC (Cascade Monte Carlo) was run 
on the neutrino-generator output. 

Above an energy threshold of 1 TeV, CMC splits cascades into several smaller, 
lower energy sub-cascades to simulate the longitudinal development of the shower. The 
sub-cascades are spaced three radiation lengths apart, where Xq ~ 41 cm for ice. Each 
sub-cascade points in the same direction as the original cascade and has a time delay 
appropriate for light travel in vacuum. The energies of the sub-cascades are taken from 
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the parameterization of the shower profile described in section 2.4. For hadronic cascades, 
CMC also applies an energy-dependent correction factor to account for the fact that hadronic 
cascades produce less Cherenkov light than their electromagnetic counterparts (see section 
2.5). 

5.3 Detector Simulation 

The result of CDRSIKA air shower simulation is a list of muons at the surface of 
the earth. For neutrino-generator muon neutrino events that have undergone a charged- 
current interaction, we also have a muon. In the first step of the detector simulation, these 
muons are propagated through the ice by the program mmc [96] . 

5.3.1 Muon Propagation with mmc 

mmc simulates the energy losses experienced by muons as they traverse the glacial 
ice at the South Pole. It takes into account ionization, bremsstrahlung, photo-nuclear 
interactions, and pair production. These energy losses in ice are summarized in figures 5.3 
and 5.4. The output of mmc is a list of muon track segments, along with their energies, and 
the secondary particles from the stochastic energy losses. 

5.3.2 Photon Propagation with photonics 

Next, we simulate the Cherenkov light emission and propagation from the muons 
and secondaries from mmc and from the electromagnetic and hadronic cascades from CMC. 
Direct tracking of the Cherenkov photons through the layered glacial ice from source to 
DOM is too computationally intensive for IceCube's Monte Carlo production. To speed up 
the process, a dedicated Monte Carlo package called photonics [84] is run beforehand to 
build lookup tables which are then used during the detector simulation, photonics tracks 
photons from a track or a cascade through the dust layers of the glacial ice, simulating the 
full wavelength-dependent scattering and absorption along the way. The result is a set of 
tables which give photon amplitudes and time delay distributions at a given location in the 
ice from a given source type and location. These tables have built into them assumptions 
about the refrozen ice which surrounds the IceCube DOM's as well as the glass transmission 
and PMT quantum efficiency. 
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Figure 5.3: Muon energy loss processes simulated by nunc. From [96]. 



5.3.3 PMT and DOM Simulation 

The simulation chain ends with the hardware simulation. For each DOM in the 
array, the photon tables are queried by a program called hit-maker to find the mean 
number of photons that particular DOM would receive from each source in the event (muon, 
radiative loss, electromagnetic cascade, etc.). A Poisson distribution with this mean is 
sampled to determine how many photons the DOM actually receives. Next, the time delay 
distributions from the photon tables are sampled to determine the arrival time of each 
photon. The result is a list of photons with their arrival times for each DOM in the array. 

Finally, this list of photons is propagated through a PMT simulation which pro- 
duces the output PMT current. This is then run through a simulation of the digitization 
electronics, local coincidence signaling, and event triggering. 

The output of an IceSim background CDRSIKA or signal neutrino simulation is a 
file identical to the experimental data produced by the detector at the South Pole. This 
file can then be run through the same processing and analysis scripts that are run over the 
experimental data. 
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Figure 5.4: Total muon energy loss as a function of energy simulated by mmc. From [96]. 
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Chapter 6 

Data Analysis 

6.1 Strategy 

In its 22-string configuration (IC-22), IceCube expected ~ 10 5 atmospheric neutrino- 
induced cascade events in its one year dataset. These signal events are buried in ~ 10 10 
downgoing muon events from cosmic ray air showers. The challenge is to develop a series of 
cuts to eliminate the background muon events while retaining sensitivity to the neutrino- 
induced cascade signal. While most of these muons should be rather easy to eliminate, we 
will eventually reach the truly difficult background for neutrino-induced cascade searches — 
muons that suffer a large radiative energy loss inside of the detector. These large radiative 
losses are indistinguishable from cascades unless there is additional early light from the 
muon track. In this chapter, we will discuss the cuts and the analysis strategies developed 
to recognize and eliminate these muons. 

Data analysis in IceCube is conducted in a blind fashion. That is, cuts are devel- 
oped using Monte Carlo simulations of signal and background, and a representative 10% 
data sample is used for validation of the background Monte Carlo. Once a set of cuts is 
developed that seems to reach the target sensitivity, the full data sample is "unblinded" 
and examined. In this way, we hope to avoid hidden biases on the part of the analyzer. 
Below we present the search strategy for atmospheric neutrino-induced cascades as it was 
developed on Monte Carlo and the 10% data sample. The next chapter contains the results 
of the full, unblinded one year dataset. 

Figure 6.1 shows a schematic of the different data processing stages, and table 6.1 
summarizes the various cuts. They will be useful guides through the rest of the chapter. 
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Figure 6.1: Schematic of the stages of the atmospheric neutrino-induced cascade analysis. 
Stages that take place at the South Pole are shaded blue and those that take place in the 
Northern Hemisphere are shaded green. Also indicated is the point where this analysis splits 
off from a parallel search for a diffuse flux of high energy extraterrestrial neutrino-induced 
cascades. 



Table 6.1: Summary table of the various cuts made in this analysis. 







Rate After Cut (Hz) 


Cut Stage 


Cut 


Background 


Signal (ve + v^) 


Trigger 


SMT 8 in 5 /xs 


558.4 


(1.6 + 87.) x 10~ 4 


Online Filter 


EvalRatio>0. 109+LFVel<0.25 


14.2 


(6.0 + 38.) x 1(T 5 


Level3 


TrackFitZenith>1.27+LlhRatio<-16.2 


3.1 


(4.9 + 5.6) x 1(T 5 


Level4a 


E > 2 TeV+ParallelogramDist < 1 . 4+L3aMLP > . 4 


1.8 x i(r 3 


(1.1 + 1.4) x 10~ 6 


Final Level 1 


E>5 TeV+MLPProd>0.73+FullAllHitSPEZ<440+FullAllHitSPEReducedLlh<10 




(1.2 + 5.2) x 1(T 7 


Final Level 2 


E>10 TeV+MLPProd>0.72+FullAllHitSPEZ<440+FullAllHitSPEReducedLlh<10 




(0.5 + 2.4) x 10~ 7 


Final Level 3 


E>15 TeV+MLPProd>0.69+FullAllHitSPEZ<440+FullAllHitSPEReducedLlh<10 




(0.3 + 1.4) x 10~ 7 



GO 
cr. 
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6.2 Online Filter at the South Pole 

The IC-22 simple multiplicity trigger (SMT) required 8 hit DOM's in a time win- 
dow of 5 fj,s. Throughout the year, the SMT rate varied between 500 and 600 Hz as 
atmospheric conditions changed with the seasons and short term weather events [97]. All 
triggered data was written to tape and retrieved from the South Pole at the end of the 
Austral winter. 

In order to proceed with analysis in a more reasonable time frame, promising 
candidate events were selected by several online filters for immediate satellite transmission 
to the Northern hemisphere. These online filters were required to run quickly so that they 
could keep up with the realtime data rate. In addition, data volume had to be reduced 
significantly in order to meet satellite bandwidth quotas. 

To select a stream of interesting events for cascade analysis, I developed and im- 
plemented a general purpose online cascade filter suitable for both low energy and high 
energy cascade studies. It was based on two quickly-calculated analytic variables: the 
tensor-of-inertia eigenvalue ratio and the linefit velocity. 

6.2.1 Tensor-of-Inertia 

The tensor-of-inertia calculation makes an analogy between photoelectrons ob- 
served in DOM's and masses distributed in space [29]. First, we calculate the center-of- 
gravity (COG) of the hits: 

COG = £>*r • n 

i=\ 

where a is the number of photoelectrons recorded by a DOM at position fl , w is an arbitrary 
amplitude weighting factor, and N c h is the number of hit DOM's. 

The COG is used as the position origin for the tensor-of-inertia, which is defined 

by 

1=1 

Continuing the analogy with a rigid body, we calculate the principal moments or eigenvalues 
of the tensor-of-inertia. 
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Figure 6.2: Online filter variables, absolutely normalized. EvalRatio measures how spherical 
an event is, and LFVel measures how fast a track moves through the detector. Cuts were 
placed to keep events with EvalRatio 0.109 and LFVel<0.25. 

Geometrically, we know that cascades look spherical in a sparsely instrumented 
detector like IceCube (see figure 2.3). So, if we calculate the three eigenvalues of the tensor- 
of-inertia they should have roughly equal values. That is, the eigenvalue ratio, defined as 
the minimum eigenvalue divided by the sum of the three eigenvalues, should be large and 
close to 1/3. For a long muon track, on the other hand, we expect the rigid body to have 
one long principal axis and two very short principal axes. The minimum eigenvalue, and 
hence the eigenvalue ratio, should be very small. The distribution of this variable for 8 
hours of trigger level data is shown in figure 6.2. A cut was placed that passes events with 
EvalRatio>0.109. 

6.2.2 Linefit 

The linefit algorithm is a simple muon track-fitting algorithm. It ignores the 
geometry of Cherenkov emission and assumes a planar wavefront traveling at a velocity v 
through the detector. Essentially, a straight line is drawn through the hit DOM's: 



fi « f+v- U 

where f\ is the location of the DOM which is hit at time ij. We want to find the position r 
and velocity v of the track that minimize the chi-squared: 
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Jv ch 

i=l 

The analytical solution to this equation is [29] 



r = (fi)-v- (U) 

?= (fj ■ U) - jn) • (U) 



where angle brackets denote an average over all hits (photoelectrons), e.g. 



*hit 



Since a muon track travels at relativistic velocities through the ice, the speed of 
the linefit should be close to the speed of light. However, for a cascade, where the light 
slowly diffuses out from a single point, the linefit velocity is expected to be small. The 
distribution of this variable is shown in figure 6.2. A cut was placed that passes events with 
LFVel<0.25. 

Together, these two cuts rejected 97% of the background while passing 37% of the 
atmospheric v e events and 4% of the atmospheric Un events. The Un efficiency appears low 
because, at this level, we're including all events, even long throughgoing tracks where 
the hadronic cascade at the interaction vertex is not inside the detector. The data rate 
from the online cascade filter varied between 17 and 22 Hz throughout the year. 

As a technical note, these reconstructions were run at the South Pole on hits 
that were required to have local coincidence span one. That is, a DOM was used in the 
calculations only if one its two nearest neighbors on a string was also hit. In contrast, the 
SMT imposed local coincidence span two, reading out a DOM only if one of its two nearest 
or two next-to- nearest neighbors on a string was also hit. Tightening the local coincidence 
for the online cascade filter was found to improve the signal-to-noise ratio. 

All subsequent stages of this analysis used only events that passed the general 
purpose online cascade filter. In the standard IceCube data files, the general purpose 
cascade filter is denoted by the key "CascadeFilter" in the FilterMask. 
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6.3 Level 2 Processing in the Northern Hemisphere 

After satellite transmission to the Northern hemisphere, more time-intensive likeli- 
hood reconstructions were run on each event. This stage of processing is known as the level 
2 standard processing and was done on a collaboration- wide basis over all data streams. 

For the cascade stream, I selected several cascade-specific reconstructions that 
would allow us to further pare down the dataset. First, the cfirst center-of-gravity first 
guess algorithm was run and used to seed a likelihood reconstruction. This likelihood 
reconstruction assumed a point-like cascade hypothesis and reconstructed the vertex and 
time using only the leading edge timing information from the waveforms and a UPandel 
probability density function. It is denoted by the name "CascadeLlhVertexFit" in the 
standard IceCube data files. See chapter 4 for details on these first guess and likelihood 
reconstructions. 

Both of these reconstructions were inadvertently run on hits that were required to 
have local coincidence span one. This oversight was corrected at a later processing stage. 

Finally, every event was subjected to a single iteration likelihood reconstruction 
assuming a muon track hypothesis. This fit is denoted by the name "TrackLlhFit" in the 
standard IceCube data files. 

No cuts were applied during level 2 processing. The quantities calculated here were 
used during the cascade-specific level 3 processing, where cuts were placed and additional 
reconstructions were run. 

6.4 Cascade Level 3 Processing 

The cascade level 3 processing started from the collaboration-wide level 2 files. We 
cut on two variables that satisfied both the needs of this atmospheric neutrino analysis and 
a parallel analysis searching for a high energy extraterrestrial neutrino signal. Additional 
reconstructions were run on the surviving events. The resulting files are referred to as the 
cascade level 3 files. 

First, a cut was placed on the zenith angle of the muon track reconstruction. This 
cut throws out events which reconstruct as obvious downgoing muons from cosmic ray air 
showers. The zenith angle distributions are shown in figure 6.3 for Monte Carlo and 8 hours 
of data, both absolutely normalized (left) and normalized to the same area (right). A cut 
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Figure 6.3: Reconstructed track zenith angle at filter level, absolutely normalized (left) and 
normalized to the same area (right). A zenith angle of zero is perfectly downgoing. A cut 
was placed that keeps events with TrackFitZenith>1.27. 

was placed that passes events with TrackFitZenith>1.27. 

Second, a likelihood ratio was formed between the track fit and the cascade vertex 
fit run at level 2. The likelihood ratio was defined as the negative log of the track likeli- 
hood minus the negative log of the cascade likelihood. A muon track should have a small 
negative log-likelihood for the track fit and a high negative log-likelihood for the cascade 
fit. Therefore the likelihood ratio should be small. The opposite holds for cascades. The 
distributions are shown in figure 6.4, absolutely normalized (left) and normalized to the 
same area (right). The oscillations are an understood artifact. A soft cut was placed that 
passes events with LlhRatio>-16.2. 

Together, these two cuts rejected another 78% of the background while passing 
82% of the atmospheric v e events and 1% of the atmospheric events. Again, the 
efficiency appears low here because, at this level, the v„ events are still dominated by long 
throughgoing tracks where the hadronic cascade at the interaction vertex is not inside the 
detector. 

After these cuts were placed, additional reconstructions were run. The first was 
the best-performing likelihood reconstruction of the cascade vertex. This uses all of the 
photoelectrons from the full set of DOM's at local coincidence span two. It's denoted by 
the name "FullAllHitSPECascadeLlh". 

Next, a split cascade reconstruction was run. The photoelectrons were split into 
two sets based on arrival time — an early set of hits and a late set of hits. The cascade vertex 
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Figure 6.4: Likelihood ratio at filter level, absolutely normalized (left) and normalized to 
the same area (right). A cut was placed that keeps events with LlhRatio>-16.2. 

for each set was individually reconstructed. For a cascade which originates at a single point 
in the detector, the two sets of hits should form an inner sphere and an outer spherical shell 
around the vertex. Therefore both vertex reconstructions should reconstruct to the same 
point in space. 

Finally, a 32-fold iterative likelihood reconstruction with a track hypothesis was 
run. This is the best track reconstruction and is denoted by "TrackLlhFit32" . 



6.5 Atmospheric Cascade Level 3a Processing 

At this point, this atmospheric neutrino analysis and the parallel high energy 
extraterrestrial analysis split and developed separate cuts. For my atmospheric analysis, I 
ran one more additional set of reconstructions on all events to produce a set of level 3a files 
and level 3a ROOT files. 

Several important reconstructions were run at this level. The first was the paral- 
lelogram distance calculation, which is a measure of how contained a reconstructed cascade 
vertex is within the detector. This is described in section 6.5.1. Next, I ran the fill-ratio 
calculation, which is described in section 6.5.2. Finally, I ran an analytic cascade energy 
reconstruction that takes into account the depth variation of the optical properties of the 
ice ("AtmCscdEnergyReco"). This energy reconstruction was developed specifically for the 
atmospheric neutrino analysis and was described in detail in chapter 4. 

In addition, several basic quality parameters were calculated for each event. Of 
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particular importance is the number of direct hits. A direct hit is defined to be one that 
arrives in a time window of -15 ns to 250 ns around the direct travel time from the cascade 
vertex to a DOM. High quality cascades which have well-reconstructed vertices generally 
have a large number of direct hits. On the other hand, muons which have a large stochastic 
energy loss are likely to leave early hits. These early hits pull the vertex time reconstruction 
toward them. Consequently, most of the hits show up as late rather than direct. 

6.5.1 Parallelogram Distance 

In AMANDA, the radial distance of a reconstructed cascade vertex from the center 
of the detector was traditionally used as a measure of containment. However, the lack of 
cylindrical symmetry in IC-22 reduces the utility of such a measure. To better capture the 
string layout of IC-22, I defined a new quantity which I denote the "parallelogram distance" . 

Figure 6.5 illustrates the construction. It shows a top view of IC-22. The solid red 
circles are the locations of the 22 strings, and the density plot shows reconstructed cascade 
vertex positions for data. 

The construction begins by defining a parallelogram that hits the innermost layer 
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of 6 strings in IC-22. This parallelogram is shown in red in the figure. By definition vertices 
that fall on this parallelogram have parallelogram distance equal to one. 

Next, for a given vertex position, this parallelogram is stretched until it hits the 
vertex. The stretching is parameterized by a single number a. Let Y u (Yj) denote the 
y-coordinate of the midpoint of the upper (lower) line of this innermost parallelogram and 
let X[ (X r ) denote the x-coordinate of the midpoint of the left (right) line of this innermost 
parallelogram. The transformation takes 



Y u -► Y u + 2(q - 1)Y U 
Y t ^Y t + 2(a - m 
Xt -> aX t 



X T 



aX r 



The transformation is defined in such a way that integer values of the parallelogram distance 
correspond to string layers of IC-22 (in general the factor need not be an integer though). 
The blue parallelogram in the figure corresponds to parallelogram distance equal to two 
and hits the main outer string layer. Parallelogram distances three and four are in green 
and magenta in the figure and hit outlying strings. 

In general, we expect that our analysis will need to focus on well-contained cascade 
events. These would be events with parallelogram distance less than 1.5 or so, as will be 
determined at later stages of the analysis. 

The distribution of this parameter at level 3a is shown in figure 6.6, absolutely 
normalized (left) and normalized to the same area (right). 

6.5.2 Fill-Ratio 

The fill-ratio calculation is intended to exploit the spherical versus track topology 
of cascades and muons. It begins with a reconstructed cascade vertex and calculates the 
distance from this vertex to all of the hit DOM's in the event. It builds the distribution of 
these distances and calculates the mean and RMS distance from the vertex. 

Next, a sphere is drawn around the reconstructed vertex. The sphere has a radius 
which is a multiple of the mean or RMS hit distance. The ratio of the number of hit DOM's 
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within the sphere to the total number of DOM's within the sphere is calculated. The two 
parameters are called the fill-ratio from mean and the fill-ratio from RMS. 

The sensitivity of this cut parameter comes from two things. First, cascades are 
rather spherical while muons are long thin tracks. A cascade is therefore expected to fill a 
larger proportion of its sphere than a muon is. However, our main background consists of 
large stochastic energy losses from muons, which are also spherical. As the energy of the 
stochastic increases though, the probability of there being an extra, early hit from the track 
increases too. This extra hit pushes the sphere to have a larger radius, so the proportion 
of hit DOM's within the sphere should go down. 

The distribution of the fill-ratio from mean at level 4a is shown in figure 6.7, 
absolutely normalized (left) and normalized to the same area (right). The distribution of 
the fill-ratio from RMS at level 4a is shown in figure 6.8, absolutely normalized (left) and 
normalized to the same area (right). 
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Figure 6.7: The fill-ratio calculated from the mean hit distance at level 4a (defined below), 
absolutely normalized (left) and normalized to the same area (right). 




Fill-Ratio Fill-Ratio 



Figure 6.8: The fill-ratio calculated from the RMS hit distance at level 4a (defined below), 
absolutely normalized (left) and normalized to the same area (right). 
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6.6 Neural Networks 

After the level 3a processing, the calculation of all of our discriminating variables 
is complete. We want to find the optimal combination of these variables, perhaps in several 
stages, in order to reject the most muon background while retaining the largest amount of 
our atmospheric neutrino-induced cascade signal. To accomplish this task, we use artificial 
neural networks (ANN) as implemented in ROOT's TMVA package [98]. 

6.6.1 Neural Networks as Multivariate Classifiers 

Let's assume that our events are characterized by N discriminating variables and 
that we want to use these variables to characterize them as signal or background. We can 
think of each event as a point in an N-dimensional space, where the coordinate along each 
axis is the value of that particular discriminating variable. If our variables are selected well, 
signal events should cluster in a certain area of this N-dimensional space and background 
events in another. The classification problem then becomes a problem of drawing a surface 
in this N-dimensional space that optimally separates the signal region from the background 
region with the least amount of leakage from one side to the other. 

There are many different multivariate classifiers in the literature, and many or 
all of them have been applied to high energy physics analyses. Examples include likeli- 
hood methods, nearest-neighbor methods, support vector machines, and neural networks. 
ROOT's TMVA package is a flexible framework for training, testing, and evaluating a large 
number of these multivariate classifiers. In testing, its neural networks outperformed all 
other methods for this analysis. 

Essentially, a neural network is just a nonlinear function. Given the N input dis- 
criminating variables, it outputs a real number. If we denote the values of the discriminating 
variables by {x\,X2, ...xjv} then: 

fANN(xi,X2,-.XN) = 2/ANN G K 

This nonlinear function is constructed by considering a collection, or network, of 
artificial neurons. Each neuron integrates the information from a set of inputs and produces 
an output. Figure 6.9 shows a schematic illustration of neuron j of a neural network. The 
inputs {y 1 ^ 1 ,y l 2~ l , ...y 1 ' 1 } are connected to the neuron with weights {w 1 ^ 1 , w 1 ^ 1 , ...w l ~^}. 
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Figure 6.9: Schematic of a single neuron from a neural network. From [99]. 

The neuron uses its "synapse function" to integrate these weighted inputs. The synapse 
function can take several forms. In this analysis, a simple sum was used: 

n 

i=i 

Next, the neuron takes the integrated input from this synapse function and uses 
its "activation function" to produce an output value. Again, the activation function can 
take several forms. In this analysis, a hyperbolic tangent function was used: 

, e Xj — e~ x i 
x~ — > y ■ = = — 

e x j + e x i 

A neural network, then, is a collection of these single neurons. Each neuron's 
output can be connected to the input of other neurons in the network. In general, Q 
neurons can have Q 2 connections between them. However, it is simplest to organize the 
neurons by layers and to only allow connections from one layer to the next. This type of 
network is known as a "multilayer perceptron" and is the type used for this analysis. It is 
indicated schematically in figure 6.10. 

In the schematic, we have N = 4 discriminating variables {x±, X2, X3, x&}. They 
enter the network from the left and are connected through a set of weights to the neurons 
in the hidden layer. These neurons, in turn, are connected by another set of weights to the 
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Figure 6.10: General layout of a multilayer perceptron (MLP). From [99]. 

output neuron which produces the final output variable j/ann- 

All neural networks used in this analysis had two hidden layers, the first with N+l 
neurons and the second with N neurons where, again, N is the number of discriminating 
variables. 

The sets of weights connecting the layers of the neural network are chosen during 
the training process. A series of signal and background events (usually from Monte Carlo) 
are fed to the network. The weights are adjusted in an attempt to push the output variable 
2/ ANN as close as possible to 1 for signal events and to -1 for background events. 

6.6.2 Level 3a Neural Network 

A first neural network was trained on five variables. These variables are listed 
below by decreasing separation power: 
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• 32FoldZenith: The track zenith angle of the 32-fold iterative muon reconstruction. 
While downgoing muons were cut with a zenith cut during level 3 processing, this 
better, more time-intensive likelihood reconstruction recovers some of the previously 
mis-reconstructed downgoing muons. 

• FullAllHitSPEReducedLlh: The reduced likelihood parameter from the best per- 
forming cascade vertex reconstruction. A lower value of this variable indicates a better 
fit to the cascade hypothesis. 

• NIHitPrac: The number of DOM's which receive a single photoelectron divided by 
the total number of photoelectrons seen by all DOM's. Muons are expected to have 
more single hits and so a higher value of this parameter. 

• NDirE: The number of direct hits from the reconstructed cascade vertex. This should 
be large for a cascade. For a muon, early hits from the muon track skew the vertex 
time. A muon, then, should have a small number of direct hits. 

• SplitDistZ: The difference in z vertex positions for the two split cascade reconstruc- 
tions. It should be close to zero for cascade events, since the early hits and the late 
hits should both reconstruct to the same position. 



Figures 6.11-6.16 show these five input variables and the output neural net classifier variable 
called L3aMLP. Plots in the right column are normalized to the same area. 
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Figure 6.11: Reconstructed track zenith angle from the 32-fold iterative reconstruction at 
level 3a, absolutely normalized (left) and normalized to the same area (right). 
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Figure 6.12: Reduced likelihood from the FullAllHitSPEReducedLlh vertex reconstruction 
at level 3a, absolutely normalized (left) and normalized to the same area (right). 
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Figure 6.13: NIHitFrac at level 3a, absolutely normalized (left) and normalized to the same 
area (right). 




Figure 6.14: Number of direct hits NDirE at level 3a, absolutely normalized (left) and 
normalized to the same area (right). 
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Figure 6.15: Split reconstruction z vertex difference SplitDistZ at level 3a, absolutely nor- 
malized (left) and normalized to the same area (right). 




Figure 6.16: Level 3a neural net output variable L3aMLP at level 3a, absolutely normalized 
(left) and normalized to the same area (right). 
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6.6.3 Level 4a Neural Networks 

A cut was placed on the level 3 neural network requiring L3aMLP>0.4. In addition, 
soft cuts were placed on the reconstructed energy and the containment of the reconstructed 
cascade vertex requiring Energy>2 TeV and ParallelogramDist<1.4. We expect that we 
will eventually have to tighten both the energy cut and the containment cut at later stages, 
but by cutting some of these events out now, further neural network training should be 
more efficient. These three cuts bring the analysis from level 3a to level 4a. 

At level4a, two individual neural networks were trained separately. Because we 
have very limited training statistics, this should result in better performance. 

The first was trained on the four variables that are the most strongly correlated 
with energy and with each other. These variables are listed below by decreasing separation 
power: 

• NDirE: The number of direct hits from the reconstructed cascade vertex still retains 
discriminating power. Again, it should be larger for cascades than for muons. 

• NHit: The total number of observed photoelectrons in all DOM's. This is correlated 
with energy. 

• LlhRatio: The likelihood from the best cascade vertex reconstruction minus the 32- 
fold iterative muon track reconstruction. This tells us which hypothesis is a better fit 
to the observed data. 

• Energy: the reconstructed energy from the AtmCscdEnergyReco fit. We expect that 
it should be easier to reject muons with higher energy radiative losses because the 
muons are more likely to leave early light. 

Figures 6.17-6.21 show these four input variables and the output neural net classifier variable 
called L4aMLPl. Plots in the right column are normalized to the same area. 
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Figure 6.17: Number of direct hits NDirE at level 4a, absolutely normalized (left) and 
normalized to the same area (right). 




Figure 6.18: Total number of photoelectrons NHit at level 4a, absolutely normalized (left) 
and normalized to the same area (right). 
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Figure 6.19: Likelihood ratio at level 4a, absolutely normalized (left) and normalized to the 
same area (right). 
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Figure 6.20: Reconstructed energy at level 4a, absolutely normalized (left) and normalized 
to the same area (right). 




Figure 6.21: First level 4a neural net output variable L4aMLPl at level 4a, absolutely 
normalized (left) and normalized to the same area (right). 
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The second neural net was trained on seven variables. These variables are listed 
below by decreasing separation power: 

• FirstHistDist: The distance from the cascade vertex to the first hit in the event. For 
cascades, the first hit should be close to the reconstructed vertex. For muons which 
leave early hits and then have a large stochastic energy loss, this distance should be 
larger. 

• ParallelogramDist: A measure of the cascade vertex containment described in sec- 
tion 6.5.1. 

• FRFillRatioFromRMS: The fill-ratio calculated from the RMS hit distance de- 
scribed in section 6.5.2. 

• FullAllHitSPEReducedLlh: The reduced likelihood parameter from the best per- 
forming cascade vertex reconstruction still retains discriminating power. A lower value 
of this variable indicates a better fit to the cascade hypothesis. 

• SplitDistZ: The difference in z vertex positions for the two split cascade reconstruc- 
tions still retains discriminating power. It should be close to zero for cascade events. 

• 32FoldZenith: The track zenith angle of the 32-fold iterative muon reconstruction 
still retains discriminating power. 

• FRFillRatioFromMean: The fill-ratio calculated from the mean hit distance de- 
scribed in section 6.5.2. 

Figures 6.22-6.29 show these four input variables and the output neural net classifier variable 
called L4aMLP2. Plots in the right column are normalized to the same area. 
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Figure 6.22: FirstHitDist at level 4a, absolutely normalized (left) and normalized to the 
same area (right). 




Figure 6.23: ParallelogramDist at level 4a, absolutely normalized (left) and normalized to 
the same area (right). 
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Figure 6.24: FRFillRatioFromRMS at level 4a, absolutely normalized (left) and normalized 
to the same area (right). 
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Figure 6.25: FullAllHitSPEReducedLlh at level 4a, absolutely normalized (left) and nor- 
malized to the same area (right). 
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Figure 6.26: SplitDistZ at level 4a, absolutely normalized (left) and normalized to the same 
area (right). 
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Figure 6.27: 32FoldZenith at level 4a, absolutely normalized (left) and normalized to the 
same area (right). 
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Figure 6.28: FRFillRatioFromMean at level 4a, absolutely normalized (left) and normalized 
to the same area (right). 




Figure 6.29: Second level 4a neural net output variable L4aMLPl at level 4a, absolutely 
normalized (left) and normalized to the same area (right). 




Figure 6.30: L3aMLP at level 4a, absolutely normalized (left) and normalized to the same 
area (right). 



In addition, the neural network output from level 3a, on which a soft cut was 
placed, still has discriminating power. It's shown in figure 6.30, absolutely normalized (left) 
and normalized to the same area (right). 
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Figure 6.31: The product of L3aMLP, L4aMLPl, and L4aMLP2 at level 4a, absolutely 
normalized (left) and normalized to the same area (right). This will be our final multivariate 
classifier. 

6.7 Final Cuts 

We have three good discriminating neural network variables at level 4a: L3aMLP, 
L4aMLPl, and L4aMLP2. These three variables are shown above in figures 6.30, 6.21, and 
6.29. For the final cut variable, it was decided to take a simple product of these three 
variables. The resulting distributions are shown in figure 6.31, absolutely normalized (left) 
and normalized to the same area (right). When normalized to the same area, the agreement 
between data and Monte Carlo is excellent. 

In figure 6.31 the contribution from atmospheric is plotted separately. As 
expected, cascades from neutral-current v„ events look identical to cascades from v e events. 
However, charged-current i/u events are more complicated. These are so-called starting 
events, where the hadronic cascade at the Vn interaction vertex is inside the detector. Below 
cut values of or so they look more background-like than the v e or the neutral-current 
events. A visual scan of these events indicates that they are ones where the outgoing muon 
contributes light to the event. As the cut variable is tightened, the hadronic cascade at the 
interaction vertex dominates the light output and they become true cascade events. 

Two additional cuts were placed before optimization of this final variable. First, 
a soft depth cut was placed on the reconstructed cascade vertex to eliminate events at the 
very top of the detector. This cut required that FullAllHitSPEZ<440 m. Next, a cut was 
placed on the reduced likelihood parameter from the best cascade vertex reconstruction. 
This was found to eliminate background from two or more muons from coincident cosmic 
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ray air showers. It required that FullAllHitSPERcducedLlh<10. 

6.7.1 Optimization 

We want to determine the optimal cut in energy and final multivariate output 
variable. For each potential energy cut, we plot a cumulative distribution that shows the 
absolute number of surviving signal and background events expected for 275.46 days of 
livetime (the full year dataset) as a function of the final multivariate cut variable. For each 
plot, CORSIKA background is normalized up to data at a cut value of 0.4. The goal is 
to look for a transition from a background-dominated region to a signal-dominated region 
before we run out of signal. However, as we will see, the background Monte Carlo statistics 
become the limiting factor. 

For an energy cut of 5 TeV, the final level distributions are shown in figure 6.32. 
Two diagnostic plots are also shown in figure 6.33. The first is a cumulative plot of the raw 
number of weighted CORSIKA background Monte Carlo events. This tells us the number 
of events before we have applied their weights. The second is a zoomed-in histogram of the 
weighted CORSIKA trace after the weights have been applied. 
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Figure 6.33: CORSIKA background Monte Carlo raw number of events (left) and sum of 
event weights (right) for an energy cut of 5 TeV. 

The cumulative distribution in figure 6.32 indicates a promising trend — the data decreases 
towards the signal prediction as the final cut variable is tightened. However, the diagnostic 
plots tell us that there is a problem with the background Monte Carlo prediction. In the 
region of highest significance, where the data is approaching the signal, we have only 2 or 
3 surviving weighted background Monte Carlo events. This causes large fluctuations in the 
background prediction, which can be seen in the large steps in the second diagnostic plot 
and which come from the loss of individual high weight events. These large fluctuations 
mean that we won't be able to accurately predict the background from our Monte Carlo, 
and a rigorous statistical detection will not be possible. 

The same conclusions hold for energy cuts of 10 TeV and 15 TeV — the data ap- 
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Figure 6.34: Final variable distribution (left) and cumulative distribution (right) for an 
energy cut of 10 TeV. 
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Figure 6.35: Final variable distribution (left) and cumulative distribution (right) for an 
energy cut of 15 TeV. 

proaches the signal prediction as the final cut variable is tightened, but the diagnostic 
plots show that we don't have enough background Monte Carlo to accurately predict the 
background. The plots are displayed in figures 6.34-6.37. 
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Figure 6.36: CORSIKA background Monte Carlo raw number of events (left) and sum of 
event weights (right) for an energy cut of 5 TeV. 
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Figure 6.37: CORSIKA background Monte Carlo raw number of events (left) and sum of 
event weights (right) for an energy cut of 5 TeV. 
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6.7.2 Final Analysis Plan 

With such limited statistics, the CORSIKA background Monte Carlo can really 
only be treated as a guide. The 10% validation data sample is the better background 
predictor. The signal prediction, coming from a known calibration source, is relatively 
robust. Looking at the data sample and the signal simulation, the trend is clear: the data 
should reach a signal dominated region as the final multivariate cut variable is tightened. 
Accordingly, we proposed to unblind three sets of cuts: 

• Final Level 1: E>5TeV + MLPProd>0.73 

• Final Level 2: E>10TeV + MLPProd>0.72 

• Final Level 3: E>15TeV + MLPProd>0.69 

After unblinding the full year's dataset with these cuts, we planned to examine the full 
cumulative distributions to look for an obvious transition from background-dominated be- 
havior to signal-dominated behavior beyond the final cut points. 

Permission was granted from the full IceCube collaboration to proceed with this 
analysis plan. The next chapter presents the results from the full 2007-2008 IC-22 dataset. 
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Chapter 7 

Results 

After a lengthy discussion and review by the IceCube collaboration, approval was 
granted to proceed with the analysis proposal outlined in the last chapter. Because back- 
ground Monte Carlo statistics were severely limited, and because the enormous computa- 
tional resources to produce more on a reasonable timescale did not exist, we decided to 
compare the observed data to the atmospheric neutrino-induced cascade signal prediction 
to look for a transition from a background-dominated region to a region where the data 
qualitatively agreed with the signal prediction. 

7.1 Bartol Atmospheric Flux Comparison 

Figure 7.1 shows the unblinded distributions and cumulative distributions of the 
final multivariate cut parameter for reconstructed energy cuts of 5 TeV, 10 TeV, and 15 
TeV. These plots make it clear that the full one year dataset does not transition to the signal 
prediction around the final cut values anticipated in the last chapter. However, zooming 
into the distributions we can see that the data does indeed start to approach the cascade 
signal prediction around a final multivariate classifier value of 0.9. Figure 7.2 shows these 
zoomed-in cumulative distributions. Steps in the distributions indicate the loss of a single 
event as the cut is tightened. 
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Figure 7.1: Final multivariate classifier distribution (left) and cumulative distribution show- 
ing the number of events surviving beyond a given classifier cut (right) for the full IC-22 
dataset with a 5 TeV (top), 10 TeV (middle), and 15 TeV (bottom) energy cut. 



122 




10 2 



10 - 



101 



0.8 



- I 1 1 1 1 1 1— 



Total Atm.v a +v^ 
Data (100%) 




0.85 



0.9 



0.95 



1.05 

Cut 



Z io 2 



■ ii ii i i i 



101 



0.8 



0.85 



0.9 



Total Atm.v B +v ( , 
Data (100%) 




0.95 



1.05 



Cut 



Figure 7.2: Zoomed-in cumulative distribution showing the number of events surviving 
beyond a given classifier cut for the full IC-22 dataset with a 5 TeV (top), 10 TeV (middle), 
and 15 TeV (bottom) energy cut. 
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Figure 7.3: The red band shows the best-fit to the AMANDA atmospheric neutrino flux 
measurement. In blue, labeled "Barr et al.", is the canonical Bartol atmospheric neutrino 
flux. The blue band is a flux measurement by Super-Kamiokande. From [101]. 

7.2 AMANDA Flux Comparison 

In comparing the observed data to the atmospheric neutrino-induced cascade sig- 
nal prediction, we need to take into account the uncertainties in the atmospheric neutrino 
flux. Up to this point, all plots in this dissertation have used the canonical Bartol atmo- 
spheric neutrino flux [19, 20] for signal normalization. In [100], the uncertainties in this flux 
were evaluated up to energies of 1 TeV. At 1 TeV, the uncertainty, which comes mainly 
from knowledge of the primaries and of the hadronic interactions, is 40% at la. Recently, 
IceCube's predecessor array AMANDA actually measured the atmospheric neutrino flux at 
the energies relevant for this analysis and found it to be higher than the Bartol flux [101]. 
The AMANDA measurements were best fit by a flux given by 

best-fit _ -iin-i\( E \ ° ^Bartol 

d~E ~ [ ' ' ' V640 GeVy dE 
This flux and its 90% error region are shown in figure 7.3. 

Figure 7.4 shows the zoomed-in cumulative distributions for the full one-year 
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dataset compared to the signal prediction using the best-fit AMANDA atmospheric neu- 
trino flux measurement. The 90% confidence region is approximated by a ±20% error on 
the best-fit flux. The ratio of the data to the signal prediction is also displayed. From these 
plots we can conclude that the data is much more consistent with the measured AMANDA 
atmospheric neutrino flux than with the canonical Bartol flux. The data-to-signal ratios at 
high cut value are generally 2 or better, indicating a signal to background ratio ~ 1. 
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Cut 

Figure 7.4: Zoomed-in cumulative distribution (left) showing the number of events surviving 
beyond a given classifier cut for the full IC-22 dataset and the signal expectation from the 
measured AMANDA atmospheric neutrino flux with a 5 TeV (top), 10 TeV (middle), and 
15 TeV (bottom) energy cut with Also displayed is the ratio of data to the neutrino signal 
prediction (right). 
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7.3 Summary Tables 

Tables 7.1-7.3 summarize the observed and predicted event numbers by class for 
the Bartol and AMANDA fluxes. The Bartol flux predictions have at least a ±40% error 
at la and the AMANDA flux predictions have a ±20% error at 90% confidence which are 
not displayed in the tables. 



Table 7.1: Summary table of observed and expected event numbers above 5 TeV. 







Bartol Flux 


AMANDA Flux 


Cut 


Observed 


v„ CC 


v„ NC 




Total v 


^ cc 


v„ NC 




Total v 


0.90 


12 


2.87 


1.61 


1.15 


5.63 


3.82 


2.15 


1.50 


7.48 


0.91 


11 


2.61 


1.45 


1.05 


5.11 


3.48 


1.94 


1.37 


6.80 


0.92 


8 


2.26 


1.27 


0.95 


4.49 


3.02 


1.71 


1.25 


5.98 


0.93 


7 


2.06 


1.14 


0.85 


4.05 


2.75 


1.53 


1.12 


5.41 


0.94 


7 


1.82 


0.96 


0.77 


3.55 


2.44 


1.29 


1.01 


4.74 


0.95 


7 


1.60 


0.85 


0.68 


3.13 


2.15 


1.16 


0.89 


4.20 


0.96 


6 


1.38 


0.74 


0.59 


2.71 


1.86 


1.00 


0.79 


3.65 


0.97 


5 


1.13 


0.66 


0.52 


2.31 


1.53 


0.90 


0.68 


3.11 


0.98 


5 


0.92 


0.56 


0.44 


1.92 


1.25 


0.77 


0.59 


2.60 


0.99 


3 


0.77 


0.52 


0.38 


1.67 


1.05 


0.71 


0.50 


2.26 


1.0 


1 


0.67 


0.43 


0.32 


1.42 


0.90 


0.59 


0.43 


1.93 
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Table 7.2: Summary table of observed and expected event numbers above 10 TeV. 







Bartol Flux 


AMANDA Flux 


Cut 


Observed 


^ cc 


^ NC 




Total v 


v„ CC 


Vp. NC 




Total v 


0.90 


8 


1.52 


0.92 


0.63 


3.07 


2.07 


1.25 


0.85 


4.16 


0.91 


7 


1.40 


0.90 


0.60 


2.90 


1.91 


1.22 


0.80 


3.93 


0.92 


6 


1.28 


0.82 


0.56 


2.65 


1.74 


1.11 


0.74 


3.60 


0.93 


5 


1.21 


0.73 


0.51 


2.45 


1.65 


1.00 


0.69 


3.34 


0.94 


5 


1.07 


0.61 


0.47 


2.16 


1.47 


0.84 


0.63 


2.95 


0.95 


5 


0.99 


0.58 


0.43 


2.00 


1.36 


0.79 


0.58 


2.73 


0.96 


4 


0.84 


0.50 


0.39 


1.73 


1.15 


0.69 


0.52 


2.37 


0.97 


3 


0.70 


0.47 


0.35 


1.52 


0.96 


0.65 


0.47 


2.09 


0.98 


3 


0.61 


0.42 


0.31 


1.34 


0.84 


0.58 


0.42 


1.84 


0.99 


2 


0.52 


0.40 


0.27 


1.19 


0.72 


0.55 


0.37 


1.63 


1.0 


1 


0.45 


0.36 


0.24 


1.05 


0.63 


0.50 


0.32 


1.45 
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Table 7.3: Summary table of observed and expected event numbers above 15 TeV. 







Bartol Flux 


AMANDA Flux 


Cut 


Observed 


^ cc 


^ NC 




Total v 


v„ CC 


y„ NC 




Total v 


0.90 


6 


0.96 


0.62 


0.39 


1.97 


1.33 


0.85 


0.53 


2.71 


0.91 


5 


0.88 


0.60 


0.37 


1.85 


1.22 


0.83 


0.50 


2.55 


0.92 


5 


0.84 


0.57 


0.34 


1.76 


1.16 


0.79 


0.47 


2.42 


0.93 


4 


0.79 


0.55 


0.32 


1.66 


1.09 


0.75 


0.44 


2.28 


0.94 


4 


0.74 


0.45 


0.30 


1.49 


1.02 


0.63 


0.41 


2.06 


0.95 


4 


0.68 


0.43 


0.28 


1.39 


0.95 


0.60 


0.38 


1.92 


0.96 


3 


0.59 


0.38 


0.26 


1.22 


0.82 


0.52 


0.35 


1.70 


0.97 


2 


0.50 


0.35 


0.23 


1.08 


0.70 


0.49 


0.32 


1.50 


0.98 


2 


0.43 


0.31 


0.21 


0.95 


0.60 


0.43 


0.29 


1.32 


0.99 


1 


0.36 


0.30 


0.18 


0.84 


0.50 


0.41 


0.25 


1.17 


1.0 


1 


0.32 


0.26 


0.16 


0.74 


0.45 


0.37 


0.22 


1.03 
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We can get a rough, intuitive feel for the relative ratios of the different neutrino- 
induced cascade signal classes displayed in the previous tables with some simple arguments. 
To have the same visible energy as a given electron neutrino charged-current event, a muon 
neutrino neutral-current event must come from a higher energy neutrino because of the 
hadronic visible energy correction (see section 2.5). This correction is around 0.8 at 1 TeV, 
so the muon neutrino must be a factor 1.25 higher in energy. The flux falls as E~ 5 ' 7 , 
but the muon neutrino flux is a factor ~ 20 higher in normalization at these energies. 
In addition, the neutral-current cross section is 1/3 of the charged-current cross section. 
Putting all these factors together, we expect roughly 2 times more muon neutrino neutral- 
current events than electron neutrino charged-current events. The ratio from the full Monte 
Carlo is ~1.4, which confirms this. 

In the same fiducial volume, muon neutrino charged-current events should occur 
3 times more frequently than muon neutrino neutral-current events because of the ratio of 
cross sections. However, this is an upper bound, since many of these charged-current events 
will have bright outgoing muons and so will not survive our cascade cuts. The ratio from 
the full Monte Carlo is around 1.8, which again makes sense given our rough estimate. 

7.4 Event Images and Waveforms 

From the cumulative distributions and tables, we see that the number of observed 
events is slightly higher than, but relatively close to, the neutrino signal prediction. Next, 
we would like to examine the 12 surviving events above the final cut value of 0.90 to see 
if they show any characteristics of muon background events. Figures 7.5-7.16 on the next 
pages show IceCube event-viewer images of these 12 surviving candidate events along with 
their reconstructed energies, arranged in order of increasing "quality" as defined by the final 
multivariate cut variable. Balloon sizes are proportional to the amount of observed light, 
and the balloon color indicates hit timing — from the earliest hits in red to the latest hits in 
blue. 

A detailed examination of these events did not uncover any muon-like character- 
istics. None of the events had early or out-of-time hits that were consistent with a muon 
track plus a large radiative energy loss. This was a necessary but not sufficient check to 
prove these events are indeed neutrino-induced cascades. 
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Figure 7.6: Run 109815, Event 14726645, E=11.2 TeV 
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Figure 7.8: Run 108479, Event 5613058, E=5.3 TeV 
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Figure 7.10: Run 109682, Event 6298338, E=133.7 TeV 
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Figure 7.12: Run 108888, Event 5189064, E=6.5 TeV 
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Figure 7.13: Run 110772, Event 16481119, E=31.3 TeV 
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Figure 7.14: Run 109366, Event 9638049, E=10.5 TeV 
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Figure 7.16: Run 109183, Event 5286833, E=15.6 TeV 
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Figure 7.17: Run 109682, Event 6298338, E=133.7 TeV with normal scaling (left) and 
shrunken scaling to show details of the innermost DOM's (right). The balloon sizes are 
proportional to the amount of detected light. 

Of these 12 events, one is a very intriguing outlier. This event is Run 109682 
Event 6298338, which reconstructs at an energy E=133.7 TeV. The event is reproduced 
in figure 7.17 along with an image where the balloon sizes, which are proportional to the 
amount of detected light, have been scaled down. The down-scaled image illustrates that 
several DOM's very close to the vertex received a tremendous amount of light in the event 
(these DOM's could not be resolved in the first image). 

Many different cross-checks were performed on this event to make sure that it is 
indeed a genuine physics event and not a detector artifact. First, monitoring and detector 
records showed that the detector was in a stable operating mode for this run — nothing out 
of the ordinary was reported and all diagnostic distributions look normal. No in-situ light 
runs were scheduled anywhere near this run, either before or after. 

To exclude the possibility of an accidental discharge of an in-situ light source, a 
dedicated "flasher" run was taken where we flashed all 12 in-situ LED's in DOM 39-19, 
which is the closest DOM to the reconstructed cascade vertex for this event. Figure 7.18 
compares the 133 TeV event to one of these flasher events. 

Several features distinguish the flashers from the 133 TeV event. First, there is a 
pronounced up-down asymmetry in the flasher events. The DOM's above the flasher see 
much more charge than do the DOM's below. For the 133 TeV event, however, DOM's 
39-19 and 39-20 see almost the same amount of charge. 
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Figure 7.18: Run 109682, Event 6298338, E=133.7 TeV compared to an in-situ flasher 
event. 

Next, the pulses in the 133 TeV event are much narrower in time than even the 
narrowest flasher pulse. These waveforms are reproduced in figure 7.19 for a flasher event 
with the broadest and the narrowest current pulse to the LED's. 

Finally, the timing is not consistent with an in-situ light source. For the flashers, 
there is a pattern to the hits: first the DOM above the flasher is hit, then the one below the 
flasher. Next comes the DOM two places above the flasher and then the DOM two places 
below the flasher. If two DOM's are neighbors, the time difference between them is never 
less than the direct travel time for 17 m. This, of course, is forbidden by causality. 

However, for the 133 TeV event, DOM 39-20 is hit first, and its leading edge time 
is at 10,000 ns plus ~ 5 ATWD samples. DOM 39-19's leading edge time is at 10,014.2 
ns plus ~ 5 ATWD samples. The time difference of ~ 14 ns is far below the direct travel 
time of 75 ns for 17 m. This means that the source of the light could not have been at 
the position of a DOM. It had to have either been between DOM's 39-19 and 39-20 on the 
string, or else somewhere off of the string entirely. So, we can conclude that the 133 TeV 
event could not have been a flasher, even an accidentally discharged flasher. 
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Figure 7.19: ATWD waveforms for flasher events compared to the 133 TeV data event. A 
flasher event at the widest setting (left), a flasher event at the narrowest setting (center), 
and the 133 TeV data event (right) are compared for DOM's 39-18 (top), 39-20 (middle), 
and 39-21 (bottom). The waveform bin width is 3.3 ns in all plots. 
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These various lines of evidence indicate that this 133 TeV event is truly an in-ice 
particle physics event and not a detector artifact or an accidental discharge of light from a 
DOM or an in-situ light source. 

All experience from this analysis as well as previous cascade analyses suggests that 
it is highly unlikely for a muon with energy greater than 133 TeV to have a large radiative 
energy loss inside the detector without first leaving early hits on an outside string. However, 
until further background Monte Carlo generation is completed, this statement cannot be 
made more quantitative. 

An obvious question that arises is the following: How likely is such an event from 
the known atmospheric neutrino flux? The expected number of events above 133 TeV for 
the livetime of this dataset and for a final cut greater than 0.9 is 0.02f At + 0.006f e = 0.03. 
The probability of a Poissonian with a mean of 0.03 fluctuating up to 1 or more counts is 
2.9%. So this event is not incompatible with the conventional atmospheric neutrino flux. 

7.4.1 Double Pulse Waveforms: A Tau Neutrino? 

Finally, there is one additional, tantalizing feature of this outlier event. The wave- 
forms from the two DOM's which receive the most charge in the event are displayed in 
figure 7.20. Both waveforms have a shoulder, a hint of a second pulse after the main pulse. 
This "double pulse" signature is one of the main signatures of a tau neutrino double bang. 
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Figure 7.20: ATWD waveforms for the 133 TeV event for the first, main DOM launch (left) 
and the second, later DOM launch (right) for DOM's 39-19 (top) and 39-20 (bottom). The 
waveform bin width is 3.3 ns. 
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This shoulder is not present in any other DOM's in the event. It's also not present 
in the flasher waveforms, which suggests that it's not an artifact of ice layering. The two 
pulses are separated by approximately 50 ns. To check the feasibility of a tau neutrino 
double bang, we can estimate the energetics of the tau track. In the best case scenario (the 
one with the lowest energy tau), the tau track points straight towards the DOM and the 
two bangs are separated by 50 ns. To live this long before decaying, the tau would have 
needed a gamma factor of 7 = St/r = 50/(2.9 x 10~ 4 ) ps 1.7 x 10 5 . This leads to a tau of 
energy E = ^m T = (1.7 x 10 5 ) • 1.776 GeV « 300 TeV. For an E~ 2 spectrum above 100 
TeV with our final cuts, our energy reconstruction has a resolution of 0.3 in log(E) with an 
offset of the reconstructed energy of -0.2. The offset would transform a 133 TeV event into 
a 210 TeV event, and the resolution means the energy would lie between 149 and 297 TeV 
at one standard deviation. The energetics, then, are at least on the right scale. 

As exciting as this event is, a full assessment of the backgrounds must be con- 
ducted. It is currently unknown whether two radiative losses along a cosmic ray muon 
track can mimic a tau, or whether the longitudinal development of a high energy cas- 
cade can do so. Another possible background is a starting muon neutrino event where the 
outgoing muon immediately suffers a large radiative loss. Within IceCube, the likelihood 
reconstruction tools for tau neutrino analysis are just now being completed and tested. A 
dedicated tau analysis on the final event sample from this analysis, taking into consider- 
ation the possible backgrounds, will be conducted in the coming months by a group at 
Pennsylvania State University. 

7.5 Comparison to Other Fluxes 

Figure 7.21 compares the observed data to several neutrino flux predictions. Red 
traces show the conventional Bartol atmospheric flux. Cyan and magenta traces show the 
sum of the Bartol conventional flux plus the contribution from two prompt neutrino models 
[21, 22]. Prompt neutrinos do not have a resolvable effect within the statistics of this 
dataset, but should become resolvable as IceCube increases its instrumented volume. 

Finally, the yellow trace shows the prediction for an extraterrestrial E~ 2 neutrino 
flux at the level of the current best neutrino-induced cascade limit from 5 years of AMANDA 
data [102]. The trace is the sum of all muon, electron, and tau neutrino contributions to 
the cascade signal. It is clear that our data excludes an extraterrestrial flux at this level 
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and that this analysis is at least as sensitive as the most-sensitive extraterrestrial analysis 
to date. However, limit setting is beyond the scope of this dissertation and will be handled 
by a separate, dedicated extraterrestrial neutrino-induced cascade analysis. 
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Figure 7.21: Cumulative distributions showing the number of events surviving beyond a 
given classifier cut for the full IC-22 dataset for 5 TeV (top), 10 TeV (middle), and 15 TeV 
(bottom) energy cuts. The signal traces show conventional atmospheric neutrinos (red), 
conventional plus prompt atmospheric neutrinos according to the models of [22] (cyan) and 
[21] (magenta), and for conventional atmospheric neutrinos plus an extraterrestrial E~ 2 
component at the current best limit [102] (yellow). 
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7.6 Systematic Uncertainties 

One remaining task is to quantify the systematic uncertainties on our atmospheric 
neutrino-induced cascade signal prediction. The dominant source of uncertainty is the 
absolute flux normalization discussed above (±40% on the Bartol model and ±20% on the 
measured AMANDA flux). Additional sources of uncertainty include: 

• Ice Properties: There are uncertainties in our knowledge of the scattering and 
absorption coefficients of the glacial ice as a function of depth. Past analyses have 
found that this contributes an overall uncertainty of ~ 15%. To quantify the effect 
of these uncertainties on this analysis, dedicated photon tables and cascade signal 
simulations will be generated and propagated through the analysis chain. 

• DOM Efficiencies: Our understanding of the glass, gel, and PMT efficiencies is 
based on lab measurements of a handful of DOM's. The simulation assumes an average 
value, and variations around this average are expected to be around ±10%. Again, 
to quantify the effect of this uncertainty, cascade signal simulations will be generated 
with varied DOM efficiencies and propagated through the analysis chain. 

• Cross Sections: The uncertainty on the neutrino-nucleon cross sections has been 
estimated to be ~ 5% [103]. 

• Absolute Energy Scale: The performance of the energy reconstruction was vali- 
dated using studies of in-situ flasher light sources (see section 4.4.3). However, addi- 
tional work must be conducted to show that the energy reconstruction performs the 
same way on data and simulation. 



145 



Chapter 8 

Conclusions and Implications 

With this analysis, we now have the first evidence that IceCube has the ability 
to see atmospheric neutrino-induced cascades. Above a reconstructed energy of 5 TeV, 
12 events were observed in the full dataset. The signal expectation from the canonical 
Bartol atmospheric neutrino flux model is 5.63 ± 2.25 events, while the expectation from 
the atmospheric neutrino flux as measured by IceCube's predecessor array AMANDA is 
7.48 ± 1.50 events. Quoted errors include the uncertainty on the flux only. These results are 
consistent with the atmospheric flux prediction plus a small amount of residual background 
contamination. Additional background Monte Carlo is needed in order to assess the purity 
of the final event sample. However, such a time intensive undertaking (six months or 
more on a computer cluster consisting of hundreds of cores) could not be completed on the 
timescale of this dissertation. 

The techniques developed in this analysis will prove useful for searches for astro- 
physical neutrinos that interact through the cascade channel. Already, the online filter and 
early cuts developed for this analysis have been used for a search for cascades from a high 
energy diffuse flux of astrophysical neutrinos [104]. Additional cuts developed for this anal- 
ysis are currently being used to search the full sky for cascades from high energy neutrinos 
from gamma ray bursts [105]. Work is underway to loosen the final cuts from this analysis 
to provide an event sample to search for correlations with flares from active galactic nuclei 
over the full sky. Finally, the possible tau neutrino candidate event will be the subject of a 
dedicated tau analysis in the coming months. 
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8.1 Implications for Future Detectors 

These results are also encouraging for future cascade analyses with IceCube. At the 
time of this writing, the 40-string configuration of IceCube (IC-40) has already completed 
its physics run. IC-40 has twice the instrumented volume of IC-22, and the dataset has 25% 
more detector uptime. That should give a factor of 2.25 improvement even without further 
refinement of the analysis techniques. The IC-40 analysis could expect to have around 50 
neutrino candidate events in the cascade channel. 

The work in this dissertation has also convinced the IceCube collaboration of the 
urgency of producing more background Monte Carlo. That work is already underway for 
the IC-40 analysis. 

At the time of this writing, the 59-string configuration of IceCube (IC-59) is cur- 
rently taking physics data. The data acquisition system for IC-59 was upgraded and should 
result in a significant improvement for cascade analysis. The IC-22 and IC-40 physics runs 
contained a hard local coincidence condition in the triggering. That is, DOM's sent data to 
the surface only if one of their two nearest or two next-to-nearest neighbors on the string 
also received light in the event. For IC-59, the physics run is being taken with a soft local 
coincidence condition in the triggering. That is, if a DOM receives light but none of its 
neighbors do, it sends a reduced subset of the event information to the surface instead of 
the full waveforms. This subset of information includes a time stamp and the total received 
charge. 

Monte Carlo studies conducted for this analysis suggest that soft local coincidence 
effectively lowers the threshold for muons to leave early hits on outer strings of the detector. 
This should significantly improve our ability to veto the muons with large stochastic energy 
losses that are the main background for cascade searches. 

In its final 80-string configuration, IceCube will be considerably better at detecting 
cascades. In addition to its quadrupled volume, superior uptime, and soft local coincidence 
data acquisition, its shape will lead to non-linear increases in background rejection ability. 
The amount of contained volume in IC-22, IC-40, and even IC-59 is not very large. When 
IC-80 is complete, there will be a large inner region surrounded by at least two layers of 
strings. This should dramatically improve our ability to veto muons with large radiative 
losses. Using the reconstruction and analysis techniques developed in this dissertation, IC- 
80 could expect to see many hundreds of atmospheric neutrino-induced cascades per year. 
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This should allow IceCube to probe models of charm production in the atmosphere. And 
in the event that IceCube detects astrophysical neutrinos, measurements in the cascade 
channel should allow us to probe the flavor ratio, perhaps telling us something about new 
physics or the conditions in the violent astrophysical sources that produce neutrinos. 
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Appendix A 

Oscillations and the Astrophysical 
Flavor Ratio 



With some simple arguments, we can understand how neutrino oscillations take 
an astrophysical flux ratio of 



and transform it to 

(z^ : v e : dearth w 1 : 1 : 1 

First, we're talking about distances on the scales of kiloparsecs and higher (1 kpc ~ 
3 x 10 16 km) and neutrinos on the order of hundreds of GeV or higher. The measured values 
of Am 2 are on the order of 10~ 3 and 10 -5 . So the phases of the sin terms in the oscillation 
probability are 



(up : v e : v T ) 
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These phases are extremely large, so we can safely replace the sin 2 term with its average 
value, 1/2. 
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Second, we know that #13 is very small, so we can safely take it to be zero. All 
imaginary CP-violating phases in the mixing matrix are proportional to sin #13, so that 
means that we're taking the mixing matrix to be fully real. Second, we know that sin 2 2#23 = 
sin 2 20 atm « 1, which means that S23 = = 1/2. Also, tan 2 #i 2 ~ 0.47 ps 1/2, so 

1 + tan 2 6»i2 = sec 2 9 12 ^c\ 2 = 2/3 

s\ 2 = 1/3 

Now we need to calculate the terms ^.{U^UpiUajUpj) = U a iUf3iU a jUpj (since U 
is real). For a = fx, [3 = fj, we have 

= U^U^U^U^ + UfaUfaU^Up! + U^U^U^U^ 
= (C12C23S12C23) 2 + (S23S12C23) 2 + (S23C12C23) 2 
= (C12C23S12C23) 2 + (C23S23) 2 

2 / 2 2 2 \ 2 

= C 23( C 23 S 12 C 12 + s 23j 

1112 1 11 

~~ 2^2 33 + T ~ 36 

Since S13 = 0, for a = /x, f3 = e which is the same asa = e,/J = /*we have 

Ufi% Uei Ufij U e j 

= UfjzUtfiUpiUei 

= -C12C23S12S12C23C12 

2 2 2 
— — C 12 C23S 12 

2 1 1 

~ ~3 2 3 
1 

~ ~9 

For the same reason, for a = e, (5 = e we have 
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We can now write down the m^ and v e oscillation probabilities: 
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So, if we start with two muon neutrinos and one electron neutrino, we'll have 
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Therefore the final flux ratio after oscillations is 1 : 1 : 1. 
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Appendix B 

Effective Livetime For Weighted 
CORSIKA Monte Carlo 

There's been some confusion and discussion within the collaboration lately about 
how exactly to define an effective livetime for weighted corsika. This note is my attempt to 
summarize my thinking on the subject. 

B.l Fundamental Definition 

For a given number of weighted corsika files, it's handy to define an effective 
livetime as a figure of merit. That way, one can say, for example, that a certain amount of 
weighted Monte Carlo corresponds to several hundred days of effective livetime at a given 
energy. 

The basic definition is this: The effective livetime is the amount of unweighted 
simulation you'd need to get the same relative error bars as the weighted simulation gives 
you. 

B.2 Simple Example 

For me, the easiest way to understand this is by considering a simple example. 
Suppose that I simulate T uw = 900s of detector livetime with unweighted corsika. This 
results in iV uw = 9 events passing a certain set of analysis cuts. This can be in a given bin 
when binned with respect to some quantity (e.g. energy) or summed over all bins to get a 
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total number of events. The relative error bar is given by 

VNZ 1 1 



Now, let's assume that we introduce some weighting scheme. This results in N w = 
9000 events passing our cuts, each with a weight Wi = 10~ 3 . Our prediction for the number 
of events stays the same, since we have 



Y^ W i = 9000 X 10 ~ 3 



The error bar is given by 



^2 wf = a/9000 x 10~ 6 « 0.1 
The relative error bar is given by 



0.01 



Ei wi 

As we can see, the weighting has reduced our relative error bar, which, after all, was the 
entire point of introducing the weighting in the first place. 

To determine the effective livetime of this weighted Monte Carlo, we ask ourselves 
the following question: How many unweighted events N e g would we need to give us the 
same relative error bar as we have with our weighted Monte Carlo? Setting the relative 
error bars equal to each other, we have 



For the example we're considering here, iV e fj ~ 10, 000. Now it's clear how to get 
the effective livetime: it's the time it would take to get N e s events out of an unweighted 
simulation which generates N uw events in T uw seconds. So we finally have the result we 
were after: 



w -(^te) (b,i) 
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For our example, T c r rj 10 6 s. 

This result scales with the number of weighted corsika files we have, as one would 
expect. If the sum of weights in one file is Y2i w i then for n files it's ~ n x Y2i w i- Likewise, 
if the sum of the squares of the weights for one file is £V wf, then for n files it's ps n x ^ u;?. 
So the livetime for n files is 



I1W 



In terms of I ceTr ay-speak, the weights lOj are given by 

weight x DiplopiaWeight 
TimeScale 

or simply by 

Wi = weight x DiplopiaWeight 
since the constant TimeScale cancels in the numerator and denominator of equation (B.l). 

B.3 Addendum 

There's actually a cancellation in equation (B.l) that results in a simpler expression 
for the effective livetime. Since the IceTray weights 

weight x DiplopiaWeight 
TimeScale 

are defined by 

weight x DiplopiaWeight N uw 

Wi = > — — - = rate m Hz = — — 

^— ' TimeScale T uw 

i i 

we can write equation (1) as 

2 



2 



EiW 



eff - /= — - [7T- =\ 



from the above definition of the IceTray weights. So we finally have an equivalent, simpler 
formulation of the effective livetime that doesn't rely on any unweighted Monte Carlo: 
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Weighted Corsika Effective Livetime Vs. Reconstructed Cascade Energy 



i i i l i i i i l i i i i l i i i i l i i i i l i i i i l i i i i l i i i 



l4aLivetime 



Entries 50 
Mean 7.013 
RMS 1.739 




log(E) [GeV] 

Figure B.l: This figure shows the effective livetime for 53,512 files from the weighted corsika 
dataset 1541 as a function of reconstructed cascade energy. The red trace is for the common 
cascade level3. The blue trace is for level4 of the atmospheric analysis. 



EX 



(B.2) 
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Figure B.2: Same as above, but plotted as a function of primary energy. 
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Appendix C 



Monte Carlo Event Weighting 



In order to normalize neutrino Monte Carlo distributions to the actual number of 
predicted events from a given neutrino flux for a given observation period, it is necessary 
to apply a weight to each Monte Carlo event. This appendix discusses that event weighting 
procedure. 

In general, the expression for number of expected events from a time-independent 
flux E) is given by 

N = T j d9.dE $(fi, E) x A cS (Q, E) 

where T is the observation time, f2 is solid angle, and A e ^ is the so-called effective area of 
the detector. Letting x = log E we have 



In E 1 

E = 10 x ' In £7 = x In 10 => x = => dx = dE => dE = E In 10 d(\ogE) 

In 10 £MnlO 

Assuming the flux is independent of the azimuth <f> and using the result above to transform 
to log E we have 



N = 2vrT / d(cos6) d(logE) E In 10 <I>(cos 9, log E) x A cff (cos 9, log E) 

The effective area ^4 e fr represents how well our detector captures neutrinos as a 
function of energy and angle and is determined from Monte Carlo simulations. We "throw" 
neutrinos of an energy E from an angle 6 at our detector, forcing them to interact, and see 
how many of them trigger the detector. The effective area is then defined as 
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A cS (cos8,logE) = ^trigg (cosg;logjg) x w 

m 'thrown 

where W is an interaction probability weight that encapsulates how likely our forced neu- 
trino interaction is to occur. We can re-write this expression as 

ai a 1 ip\ dN tTigg /d(cos6)d(logE) 

A eS (cos9,logE) = — — cosfl,log£ x W 

dN thrown /d{cos 0)d(log E) 

Substituting this expression for A e g into the integral for N above, the factor d(cos 9)d(\og E) 
in the numerator will cancel the variables of integration. We're left with an integral over 



dNfa jgg . 



N = 2vrT / dN tligg 



E In 10 x W x $(cos0,log£) 



diV throwil /d(cos 0)d(log E) 
It's clear now that the expression in square brackets is our event weight, which we apply 

to each triggered event before we sum over all triggered events to get the actual number of 

expected events. 

The factor in the denominator of the weight is determined by the parameters of 
the thrown Monte Carlo. Typically, this thrown Monte Carlo has either an E~ l or E~ 2 
energy spectrum. For a thrown Monte Carlo spectrum of E -1 we have 

d-^thrown ^ 



d(cos9)d(logE) 

A here is the normalization constant. Integrating over cos 9 from —1 to 1 and over logE 
from logii^min to log-Emax and requiring the total number of thrown events to be ATthrown 
we get an expression for A: 

^■^thrown ■''thrown 



d(cos#)d(log£') 2 logf^ 

^rnin 

One more subtlety remains. Monte Carlo generation typically throws an equal 
number of v and v together. If our flux is expressed in terms of the sum of v + u, this 
expression for A is fine. However, if our flux is expressed in terms of v and v separately, as 
is the case for atmospheric fluxes, then we need to add a factor of 2 to account for the fact 
that we've actually thrown Athrown/2 v and Athrown/2 v. In this case A is given by 

^-^thrown -^thrown 



d{cos9)d{logE) 2x2 log §- 
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For a thrown Monte Carlo spectrum of E we have 

^thrown B In 10 

d{cos6)d{logE) ~ E 

Again, integrating over cos# and \ogE we get an expression for B: 



q -^thrown Brain B-a 



2 E mSLX E m ' m 

and again, if the flux is expressed separately in terms of v and v we get an extra factor of 
two: 

-^thrown B m [ n E m&K 



2x2 E m&x E m \ n 



Our expression for ,, r 2tEm is thus 



d(cos 9)d(\og E) 

dNfbiavm B ln 10 ^thrown ^ 10 E min E n 



d(cos6)d(logE) E 2E E m&x -E n 

or 

diVthrown B 111 10 thrown In 10 E m - m E ma 



-tain 



d(cos 0)d(log E) E 2 x 2E £ max - E n 

We now have all the ingredients that go into the event weight to properly normalize event 
number expectations for different types of Monte Carlo. 
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